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ABSTRACT 

The positions and velocities of galaxies in the Local Group (LG) measure the grav¬ 
itational field within it. This is mostly due to the Milky Way (MW) and Andromeda 
(M31). We constrain their masses using distance and radial velocity (RV) measure¬ 
ments of 32 LG galaxies. To do this, we follow the trajectories of many simulated 
particles starting on a pure Hubble flow at redshift 9. For each observed galaxy, we 
obtain a trajectory which today is at the same position. Its final velocity is the model 
prediction for the velocity of that galaxy. 

Unlike previous simulations based on spherical symmetry, ours are axisym- 
metric and include gravity from Centaur us A. We find the total LG mass is 
4.33^0.32 with 0.14 ± 0.07 of this being in the MW. We approximately 

account for IC 342, M81, the Great Attractor and the Large Magellanic Cloud. 

No plausible set of initial conditions yields a good match to the RVs of our sample 
of LG galaxies. Observed RVs systematically exceed those predicted by the best¬ 
fitting ACDM model, with a typical disagreement of km/s and a maximum 

of 110 ± 13 km/s for DDO 99. Interactions between LG dwarf galaxies can’t easily 
explain this. 

One possibility is a past close fiyby of the MW and M31. This arises in some 
modified gravity theories but not in ACDM. Gravitational slingshot encounters of 
material in the LG with either of these massive fast-moving galaxies could plausibly 
explain why some non-satellite LG galaxies are moving away from us even faster than 
a pure Hubble flow. 

Key words: galaxies: groups: individual: Local Group - Galaxy: kinematics and 
dynamics - Dark Matter - methods: numerical - methods: data analysis - cosmology: 
cosmological parameters 


1 INTRODUCTION 

In a homogeneous universe, particles would follow a pure 
Hubble flow. This means their velocities would depend on 
their positions according to 

V (t) = H (t) r where H = Hubble parameter at time t (1) 

However, the Universe is inhomogeneous on small 
scales. The resulting inhomogeneous gravitational field 
causes motions to deviate from Equation 1. These devia¬ 
tions — termed ‘peculiar velocities’ — are easily discerned 
in the Local Group (LG). Thus, the observed positions and 
velocities of LG galaxies hold important information on the 
gravitational field in the LG, both now and in the past.^ 

* Email: ib45@st-andrews.ac.uk (Indranil Banik), 
hz4@st-andrews.ac.uk (Hongsheng Zhao) 

^ Due to Hubble drag (paragraph below Equation 22), peculiar 
velocities are mostly sensitive to forces acting at late times. 


Therefore, by investigating a range of physically motivated 
models for the gravitational field of the LG, we can hope to 
see which ones — if any — plausibly explain these observa¬ 
tions. This technique is known as the timing argument. 

The timing argument was first applied to the Milky Way 
(MW) and Andromeda (M31) galaxies over 50 years ago 
(Kahn V Woltjer 1959). This pioneering work attempted to 
match the present relative velocity of the MW and M31, 
assuming no other major nearby sources of gravity. As M31 
must initially have been receding from the MW but is cur¬ 
rently approaching it at ~110 km/s (Slipher 1913; Schmidt 
1958), it was clear that models with very little mass in the 
MW & M31 could not work.^ In fact, their combined mass 
had to be ~10 times the observed baryonic mass in these 
galaxies. This provided one of the earliest indications that 
most of the mass in typical disc galaxies might be dark. 


^ In the limit of no mass, M31 would be receding at ^^50 km/s. 
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This conclusion has withstood the test of time, at least 
in the context of Newtonian gravity. More recent works find 
a total LG mass of M ~ (4 — 5) x 10^^M© (Li & White 2008; 
van der Marel et ah 2012b; Partridge et al. 2013). This is 
roughly consistent with the combined dynamical masses of 
the MW and M31. For example, analysis of the giant south¬ 
ern stream around Andromeda (a tidally disrupted satellite 
galaxy) yielded 2.5 x lO^^M© (Fardal et al. 2013).^ 

Combining a wide variety of observations of our own Galaxy, 
McMillan (2011) found that 1.5 x IO^^Mq.^ How¬ 

ever, careful analysis of the Sagittarius tidal stream (New- 
berg et al. 2002; Majewski et al. 2003) found a mass of about 
half this (Gibbons et al. 2014), though this depends on the 
uncertain distance to the progenitor. 

The timing argument seems to suggest a higher mass 
than the sum of the MW & M31 dynamical masses. The 
tension would be further exacerbated if the LG mass was 
smaller in the past, forcing up the present mass inferred by 
the timing argument. This is quite likely as galaxies accrete 
mass from their surroundings. 

One possible explanation may be that, in the context 
of a cosmological simulation, the timing argument overesti¬ 
mates the LG mass (Gonzalez et al. 2014). However, this 
trend is not seen in the work of Partridge et al. (2013), 
whose timing argument calculations included the effect of 
dark energy. In any case, the tension does not appear to be 
significant. 

The present Galactocentric radial velocity (GRV) of 
Andromeda provides just one data point. Therefore, it can 
only be used to constrain one model parameter: the total 
LG mass. The mass ratio between the MW and M31 can’t 
be constrained in this way, although it is likely on other 
grounds that as M31 is larger (Bovy & Rix 

2013; Courteau et al. 2011) and rotates faster (Carignan 
et al. 2006; Kafle et al. 2012). 

More importantly, we can’t determine if the model itself 
works with just one data point. As a result, it has been 
suggested to include more distant LG galaxies in a timing 
argument analysis (Lynden-Bell 1981). Such an analysis was 
attempted a few years later (Sandage 1986). This work sug¬ 
gested that it was difficult to simultaneously explain all the 
data then available. 

The quality of observational data has improved sub¬ 
stantially since that time. More galaxies have also been dis¬ 
covered, providing additional constraints on any model of 
the LG. This is partly due to wide field surveys such as the 
Sloan Digital Sky Survey (SDSS, York et al. 2000) and the 
Pan-Andromeda Archaeological Survey (McGonnachie et al. 
2009). 

Such surveys have shown that satellite galaxies of the 
MW are preferentially located in a thin (rms thickness 
~25 kpc) co-rot at ing planar structure (Pawlowski & Kroupa 
2013). Known MW satellites were mostly discovered using 
the SDSS, which has only limited sky coverage. Even when 
this is taken into account, it is extremely unlikely that the 
MW satellite system is isotropic (Pawlowski 2016). In fact, 
this hypothesis is now ruled out at > 5cr. 

A similar pattern is also evident with the satellite galax- 


^ This is an estimate of M2oo- 
^ This is an estimate of the virial mass. 


ies of Andromeda (Ibata et al. 2013). Roughly half of its 
satellites are consistent with an isotropic distribution but 
the other half appears to form a co-rotating planar struc¬ 
ture even thinner than that around the MW. However, co¬ 
rotation can’t be definitively confirmed until proper motions 
become available. 

The observed degree of anisotropy appears very difficult 
to reconcile with a quiescent origin in a Lambda-Gold Dark 
Matter (AGDM) universe (Pawlowski et al. 2014, and refer¬ 
ences therein). This result seems to hold up with more recent 
higher resolution simulations (Gillet et al. 2015; Pawlowski 
et al. 2015). One reason is that filamentary infall is unlikely 
to work because it leads to radial orbits, inconsistent with 
observed proper motions of several MW satellites (Angus 
et al. 2011). 

This result has recently been challenged by Sawala et al. 
(2014) and Sawala et al. (2016) based on the eagle simu¬ 
lations, which include baryonic physics (Schaye et al. 2015; 
Grain et al. 2015). When comparing with the observed satel¬ 
lite systems of the MW and M31, these investigations did 
not take into account all of the available information, in 
particular the observed distances to the MW satellites. Once 
this is considered, it becomes clear that the observed dis¬ 
tribution of satellites around the MW is very anisotropic, 
making a quiescent scenario for their origin much less likely 
(Pawlowski et al. 2015). Moreover, the inclusion of baryonic 
physics had very little impact on the extent to which satellite 
systems are anisotropic. This is what one would expect given 
the large distances to the MW satellites. 

In this context, it seems surprising that a recent investi¬ 
gation found that the observed satellite systems of the MW 
and M31 are consistent with a quiescent AGDM origin at the 
5% and 9% levels, respectively (Gautun et al. 2015). How¬ 
ever, this analysis suffers from several problems, in partic¬ 
ular not considering several objects orbiting the MW (only 
its 11 classical satellites are considered). The result for the 
MW is based on assuming that | of the sky is not observ¬ 
able due to the Galactic disc. The actual obscured region 
is likely smaller, making the observed distribution of MW 
satellites harder to explain. Some of the more important 
deficiencies with this investigation have been explained by 
Lopez-Gorredoira & Kroupa (2016, last paragraph of page 
2 ). 

The MW and M31 are ~ 0.8 Mpc apart now (Mc¬ 
Gonnachie 2012) and have never interacted in AGDM (see 
Figure 1). Thus, one might expect their satellite systems 
to be almost independent in this model. Indeed, it has re¬ 
cently been demonstrated in simulations that the degree of 
anisotropy of the MW satellite system is not enhanced by 
the presence of an analogue of M31 (Pawlowski & McGaugh 
2014b). 

It must be borne in mind that all these authors focused 
on Local Group satellites merely because they happen to be 
nearby, allowing for much more accurate measurements of 
3D positions and velocities. It is very difficult to conduct 
similarly detailed investigations further away. Thus, while it 
may be dangerous to conclude too much about the Universe 
based on just ~ 50 satellite galaxies, one should at least 
concede that these are located in two essentially independent 
systems which were not selected because of their anisotropy. 

Although a quiescent origin for these highly anisotropic 
satellite systems appears unlikely, it is possible that an an- 
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cient interaction created them by forming tidal dwarf galax¬ 
ies (TDGs, Kroupa et ah 2005). After all, there are sev¬ 
eral known cases of galaxies forming from material pulled 
out of interacting progenitor galaxies (e.g. in the Antennae, 
Mirabel et ah 1992). 

Such TDGs tend to be more metal-rich than primordial 
galaxies of the same mass (e.g. Groxall et al. 2009). M31 
satellites in the planar system around it seem not to have 
different chemical abundances to M31 satellites outside this 
plane (Gollins et al. 2015). This might be a problem for 
the scenario, had it involved a recent interaction. But with 
a more ancient interaction, the problem seems to be much 
less severe (Recchi et al. 2015). Essentially, this is because 
gas in the outer parts of the MW/M31 would have been very 
metal-poor when the interaction occurred. This would lead 
to TDGs that were initially metal-poor, similar to primordial 
objects of the same age. 

TDGs should be free of dark matter as their escape 
velocity is much below the virial velocity of their progenitor 
galaxies (Barnes & Hernquist 1992; Wetzstein et al. 2007). 
Thus, a surprising aspect of LG satellite galaxies is their high 
mass-to-light (M/L) ratios (e.g. McGaugh & Milgrom 2013). 
These ratios are calculated assuming dynamical equilibrium. 
Tides from the host galaxy are probably not strong enough 
to invalidate this assumption (McGaugh & Wolf 2010). With 
dark matter unlikely to be present in these systems, the high 
inferred M/L ratios would need to be explained by modified 
gravity. 

One possibility is to use Modified Newtonian Dynam¬ 
ics (MOND, Milgrom 1983). This imposes an acceleration- 
dependent modification to the usual Poisson equation of 
Newtonian gravity (Bekenstein & Milgrom 1984). Despite 
having only one free parameter, MOND fares well at explain¬ 
ing rotation curves of disc galaxies (Famaey & McGaugh 
2012, and references therein). It also seems to work for LG 
satellites (McGaugh & Wolf 2010; McGaugh & Milgrom 
2013), although the relevant observations are challenging. 

Applying this theory to the MW and M31, Zhao et al. 
(2013) found that they would have undergone an ancient 
close flyby ~9 billion years (Gyr) ago. The thick disc of the 
MW would then be a natural outcome of this interaction. 
Indeed, recent work suggests a tidal origin for the thick disc 
(Banik 2014). Moreover, its age seems to be consistent with 
this scenario (Quillen Sz Garnett 2001). 

An ancient flyby of M31 past our Galaxy might have 
affected the rest of the LG as well. Infalling dwarf galaxies 
might have been flung out at high speeds by gravitational 
slingshot encounters with the MW/M31. Material might also 
have been tidally expelled from within them, perhaps form¬ 
ing a dwarf galaxy later on. As a result, the velocity field 
of the LG would likely have been dynamically heated. We 
hope to investigate whether there is any evidence for such a 
scenario. 

To this end, the use of more distant LG galaxies can 
be particularly useful. Within the context of AGDM, Penar- 
rubia et al. (2014) used non-satellite galaxies within ~ 3 
Mpc for a timing argument analysis. Satellite galaxies can’t 
easily be used in this way because the velocity field becomes 
complicated close to the MW or M31 (Figure 3). Intersect¬ 
ing trajectories make it difficult to predict the velocity of a 
satellite galaxy based solely on its position. 

We perform a similar analysis of the same ‘target’ galax¬ 


ies as in that work. The basic idea is the same: we construct a 
test particle trajectory that today is at the same position as a 
target galaxy.^ The final velocity relative to the MW is then 
projected onto our line of sight (Equation 50). This model- 
predicted GRV is corrected for the motion of the Sun with 
respect to the MW, yielding a heliocentric radial velocity 
(HRV) prediction which can be compared with observations. 
When proper motion measurements become available, it will 
be very interesting to compare the full 3D velocities of LG 
galaxies with our models. 

For simplicity, we assume that the only massive objects 
in the LG are the MW and M31, which we take to be on 
a radial orbit. Recent proper motion measurements of M31 
indicate only a small tangential motion relative to the MW 
(van der Marel et al. 2012a). This makes the true orbit al¬ 
most radial. 

The recent work of Salomon et al. (2016) argues for 
a high M31 proper motion (~100 km/s) based on redshift 
gradients in the M31 satellite system. This measurement 
is consistent with the more direct measurement of van der 
Marel et al. (2012a), though there is some tension. This 
might be explained by intrinsic rotation of the M31 satellite 
system. With a field of view of perhaps 5°, rotation at only a 
~10 km/s level can masquerade as a proper motion of ~100 
km/s. In fact, there is strong evidence that nearly half of 
the M31 satellites rotate coherently around it (Ibata et al. 
2013). Although Salomon et al. (2016) take this into account 
to some extent, other rotating satellite planes might also 
exist around M31. This is suggested by recent investigations 
into the kinematics of its globular cluster system (Veljanoski 
et al. 2014). Moreover, a large tangential velocity between 
the MW and M31 would show up as a dipole-like feature in 
the radial velocities of distant LG galaxies. This has been 
searched for but not found (Peharrubia et al. 2016). Thus, 
we assume the van der Marel et al. (2012a) proper motion 
measurement is more accurate, making the MW—M31 orbit 
nearly radial. 

Starting at some early initial time we evolved for¬ 
wards a large number of test particles in the gravitational 
field of the MW and M31. We took the barycentre of the 
LG dX t — ti as the centre of the expansion. The initial 
velocities followed a pure Hubble flow (Equation 37). This 
is because the Universe was nearly homogeneous at early 
times — peculiar velocities on the last scattering surface are 
only ~1 km/s (Planck Gollaboration 2015), much less than 
typical values today (~50 km/s, see Figure 11). 

As both the initial conditions and the gravitational 
field are axisymmetric, test particles move within meridional 
planes (i.e. those containing the symmetry axis). This al¬ 
lowed us to use an axisymmetric model. We briefly mention 
that the gravitational field in our model varies with time, 
because the MW and M31 move. 

A major improvement of our analysis is that the LG is 
not treated as spherically symmetric. This assumption is not 
a very good one as the targets considered are at distances 
of ~1—3 Mpc. Meanwhile, the MW—M31 separation is ~0.8 
Mpc (McGonnachie 2012). This means that the gravitational 
potential — and thus velocities — are likely to deviate sub- 


^ Our model is effectively two-dimensional, so we used a 2D ver¬ 
sion of the Newton-Raphson algorithm to achieve this. 
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stantially from spherical symmetry in the region of interest. 
However, we expect only small deviations from axisymmetry 
for reasons just stated. 

Objects outside the LG can have some influence on our 
results because they can raise tides on the LG. The most 
important perturbers were identihed by Peharrubia et al. 
(2014) as M81, IG 342 and Gentaurus A. Their properties are 
given in Table 1. We directly included the gravity of the most 
massive of these objects, Gen A (Section 2.2.2). We took 
advantage of its location on the sky being almost exactly 
opposite that of Andromeda. Due to the large distance of 
Gen A from the LG (~ 4 Mpc), its velocity is dominated 
by the Hubble expansion (Karachentsev et al. 2007). This 
makes the LG—Gen A trajectory almost radial, allowing us 
to continue using our axisymmetric model. 

Our paper is structured as follows: The governing equa¬ 
tions and methods are described in Section 2. This section 
also shows some results, to give a rough idea of what happens 
in our simulations. Gomparison of simulation outputs with 
observations is done in Section 3. The posterior probability 
density functions of all variables and pairs of variables are 
shown in Figure 7. Our results indicate that no model comes 
close to reproducing all the observations simultaneously. 

In Section 4, we discuss several shortcomings of our 
model and whether accounting for some of them might help 
to explain the observations. In Table 5, we show how Gen 
A affects our results. We also estimate carefully the effects 
of M81 & IG 342 (Table 6), the Great Attractor (Figure 16) 
and the Large Magellanic Gloud (Figure 17). These objects 
seem to little affect GRVs and often worsen the discrepancy 
with the best-fitting model. We suggest a possible explana¬ 
tion for our results in Section 4.6. Differences between our 
approach and the similar study of Peharrubia et al. (2014) 
are described in Section 4.7. Our conclusions are summarized 
in Section 5. 


2 METHOD 

The method we follow is to ensure a simulated test parti¬ 
cle ends up at the same position as each LG galaxy in our 
sample (a ‘target’). At present, only the radial velocities of 
our targets are available. Thus, the velocity of this particle 
relative to that of the MW is projected onto the direction 
towards the particle (Equation 50). This model-predicted 
GRV is then corrected for solar motion in the MW and 
compared with observations. The procedure is repeated for 
different model parameters, which are systematically varied 
across a grid. Therefore, within the priors we set (Table 2), 
all model parameter combinations were investigated. 


2.1 Equations of motion 

We begin with the metric in the weak field limit 


^2 2,2 

as — c dr 


s (x) = 


C{x) = 


sinh (x) in an open universe 

X in a flat universe 

cosh (x) in an open universe 

1 in a flat universe 


( 2 ) 

(3) 

(4) 

(5) 


Here, c is the speed of light and r is proper time. The 
scale-factor of the universe is a. The spatial part of the met¬ 
ric has been written in spherical polar co-ordinates, with 
dO. representing a change in angle. Using the co-ordinates 
x°=t, x^=x? = 0 and assuming spherical symmetry, 

we get a diagonal metric where 

floo = (l + ^) (6) 

flu = -a (l - ^) (7) 

922 = (l - ^)(x) = flii-S'^(x) (8) 

As the metric coefficients are independent of = 0, 
the geodesic equation tells us that 

3 

X2 E g2bX^ (only non-zero term is 5 = 2) (9) 

6=0 



= constant (10) 

We use q to denote the derivative of any quantity q with 
respect to proper time. In weak gravitational fields (4> ^ c^), 
proper and co-ordinate time are almost equal, making r ~ t. 
Bearing this in mind. Equation 10 tells us that the specific 
angular momentum of a test particle is conserved. This is 
due to the spherical symmetry of the situation. 

Eor the radial component of the motion, we use the 
geodesic equation in the form 
3 3 

-^+EE r^bcX^x^ = 0 where (11) 

b=0 c=0 

i^bddc ^c9bd ^ddbc) ( 12 ) 

^ d=0 

Here, we use the notation d^^q = for any quan¬ 
tity q. The non-zero Ghristoffel symbols relevant to a non- 
relativistic test particle in this situation are 


r'oo 

4)' 

(13) 

r'oi 

a 

(14) 

rSa - 

^ -^(x)C(x) 

(15) 

r\i p 

4>' 

(16) 


Here, q implies a partial derivative with respect to the 
co-moving co-ordinate x rather than physical distance ax- 
Putting in the non-negligible Ghristoffel symbols^ into Equa¬ 
tion 11, we get that 

X+% + 2Hx-S{x)C{x)e'' = 0 (17) 

^ The r^ii term effectively causes 4>' —)■ <!>' ^1 — where v 
is the speed of the particle with respect to a co-moving observer 
at the same place. This leads to a special relativistic correction 
which makes it difficult for a potential gradient to accelerate a 
particle if its speed is close to that of light. For non-relativistic 
particles, the effect of this term is negligible because v ^ c. 
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In terms of physical co-ordinates r = ax, this becomes 


r = (^^x + ‘2Hx + x'^ a (18) 

/a 1 

= ^r-^ + Six)C{x)ae" ( 20 ) 

a or 

The real Universe is close to spatially flat (Planck Col¬ 
laboration 2015). Thus, from now on, we will only consider 
the case of a flat universe. This is defined as one having a 
density p equal to the critical density pcrit, if we count both 
matter and dark energy towards the total density. 

P — Pcrit = (21) 

Equation 21 is valid at all times, although both p and 
H vary with time. In such a universe. Equation 20 becomes 


r 


a 

— r — 
a 


^ +r6> 
or 


( 22 ) 


This looks very similar to the Newtonian equation of 
motion. The last term corresponds to the centrifugal force 
while the ^ term corresponds to the potential gradient. 
The only novel aspect is the term ^r. The importance of this 
term becomes clear if we consider a homogeneous universe, 
meaning that $ = 0 everywhere and at all times. In this 
case, we expect the distance between two non-interacting 
test particles to behave as r oc a (t) (i.e. their co-moving 
distance is constant). This implies that r = ^r. This term 
has also been called Hubble drag because it tends to re¬ 
duce the magnitude of peculiar velocities.^ In the absence 
of potential gradients, we would get Vpec oc where the 
peculiar velocity is defined by 


Vpec = r — Hr 


(23) 


In general, the Universe is neither homogeneous nor 
spherically symmetric. Eor such circumstances, we suppose 
that the generalization of Equation 22 is given by 


r — —r — 
a 


(24) 


With the equations of motion in hand, we now need to 
relate the potential $ to the density perturbations that act 
as its source. To do this, we use the 00 component of the 
field equation of General Relativity. 


R., = where (25) 

3 

T = (26) 

d=0 


Here, is the energy-momentum tensor while is 
the Ricci tensor, related to the curvature of the metric. Per¬ 
turbations to the solution for a homogeneous universe must 
satisfy the equation 

^Roo = -^A (too - irg„„) (27) 


^ If the Universe were contracting, then this term would be a 
forcing to peculiar velocities rather than a drag upon them. 


In this case, for non-relativistic sources which are almost 
pressureless (like baryons and cold dark matter), we get that 

Rqo ~ ~ (sum over spatial indices only) (28) 

i=l ^ ^ 

The stress-energy tensor takes on a particularly simple 
form: its only non-zero element is Tqq pc^. Thus, 

floo^ ~ Ho ~ (29) 


Using Equations 28 and 29 in Equation 25, we get that 


= 47rGAp 


(30) 


Here, is the Laplacian of 4> with respect to physical 
co-ordinates. In spherical symmetry, it is 




! 


dr 


(31) 


Equation 30 is very similar to the usual Poisson equa¬ 
tion of Newtonian gravity. Note, however, that only devia¬ 
tions from the background density act as a source for $ (i.e. 
it is sourced by Ap rather than p). 


2.2 Simulations 


2.2.1 Including the Milky Way and Andromeda 


The LG is assumed to consist of two point masses (the MW 
and M31) plus a uniform distribution of matter at the same 
density as the cosmic mean value p (see Section 4.1 for fur¬ 
ther discussion of this point). Our simulations start when 
the scale factor of the Universe a — a^. We used a^ — 0.1, 
though our results change negligibly if a. = instead (see 
Section 4.6). 

The initial separation of the MW and M31 di is var¬ 
ied to match their presently observed separation do using 
a Newton-Raphson technique. Note that altering di alters 
their initial velocities because the galaxies are assumed to 
have zero peculiar velocity at the start of the simulation 
{t — ti). Thus, their final attained separation df depends 
strongly on di. 

The MW—M31 orbit is taken to be radial, a reason¬ 
able assumption given their small tangential motion (~17 
km/s compared to a radial velocity of ~110 km/s, van der 
Marel et al. 2012a). This makes the gravitational field in the 
LG axisymmetric. As the initial conditions are spherically 
symmetric (Equation 37), a 2D model is sufficient for this 
investigation. 

Applying Equation 24 to a radial orbit, the distance d 
between the galaxies satisfies 


d 

d 


GM a , 

'^i2~ + 
d^ a 


id. di initially 


(32) 

(33) 


H is the value of the Hubble constant - when t — ti, 
while di is the MW—M31 separation at that time. M is 
the combined mass of the MW & M31. It can be verified 
straightforwardly that when M = 0 (i.e. non-interacting test 
particles), we recover d oc a. In this case, the galaxies trace 
the cosmic expansion but don’t influence each other. 

Equation 32 implicitly assumes that the MW and M31 
are surrounded by a distribution of matter with the same 
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-6 -4 -2 0 2 4 6 8 

Fractional difference in final between iterations ^ 10'^ 


Figure 1. MW—M31 separation d(t) for a typical case where 
= 0.2 and Mi = 3.4 x 10^^Mq. d{t) always looks broadly 
similar — in ACDM, the MW and M31 have never approached 
each other closely for any plausible model parameters. 


Figure 2. Fractional difference between the final attained MW 
mass on the first and second runs of each simulation. The very 
small values show that our solution for (t) converged well. 

Simulations without Centaurus A give similar results. 


density as the cosmic mean value p. This point is discussed 
more thoroughly in Section 4.1 where we also redo our entire 
analysis assuming instead that the surroundings of the MW 
and M31 are empty. 

We use a standard flat^ dark energy-dominated cosmol¬ 
ogy with parameters given in Table 2. Therefore, 

d 47rG . ^ X 

— — - {pm — 2pa) (34) 

= + f^A.o) (35) 

Defining time t to start when a = 0 and requiring that 
d = Hq when a = 1 (the present time), we get that 

a{t) = (36) 

The present values of Hq and Qm,o uniquely determine 
the present age of the Universe tf via inversion of Equation 
36 to solve for when a = 1. We also use it to determine when 
a = , thereby fixing the start time of our simulations. 

The timing argument is particularly sensitive to late 
times (Figure 4). This makes it important to correctly ac¬ 
count for the late-time effect of dark energy. Because this 
tends to increase radial velocities of LG galaxies, one is 
forced to increase the mass of the LG to bring their pre¬ 
dicted radial velocities back down to the observed values. 
As a result, the inclusion of dark energy in timing argument 
analyses of the LG increases its inferred mass by a non- 
negligible amount (Partridge et al. 2013). 

Once we obtained a trajectory that (very nearly) satis¬ 
fied df = d{tf) = do, we had the ability to find the gravi¬ 
tational field everywhere in the Local Group at all times. A 
large number of test particles were then evolved forwards, all 
starting on a pure Hubble flow with the centre of expansion 
at the barycentre of the Local Group. 

V, = (37) 


T Ua,0 — 1 


Note that Equation 37 also applies to the MW and M31, 
which we model as point masses. A point mass approxima¬ 
tion should work for determining d(t) as Andromeda never 
gets very close to the Milky Way (Figure 1). However, it 
is not good for handling close encounters of test particles 
with either galaxy. Thus, we adjust the forces they exert on 
test particles to be (x ^ at low r (i.e. close to the attracting 
body). This is for consistency with the observed flat rotation 
curves of the MW and M31. To recover ^ (x ^ at large r, 
we set the gravity towards each galaxy to be 


9 

b 


GM b 

y/TTF 


where 


r 


s 


(38) 

(39) 


is chosen so that the force at r ^ leads to the 
correct flat line level of rotation curve for each galaxy, i.e. 

MW, we take — 180 km/s (Kafle 

et al. 2012) while for Andromeda, we take — 225 km/s 
(Garignan et al. 2006). 

Gombining Equation 38 with the cosmological acceler¬ 
ation term, the equation of motion for our test particles is 


- ^ ^ 

j=Mw, M31 [\r — r 


GMj {r-r.) 




(40) 


Some trajectories go very close to the MW or M31. Ap¬ 
proaches within a distance of race (given in Table 2) are 
handled by terminating the trajectory and assuming the 
particle was accreted by the nearby galaxy. This causes the 
mass of that galaxy to increase. 

As we solved the test particle trajectories sequentially, it 
wasn’t possible until the very end to have the mass histories 
Mmw{^) ^M 3 i(^)- Thus, we assumed constant masses 
for the force calculations. We then repeated the process, us¬ 
ing the previously stored mass histories for each galaxy. This 
meant that the initial MW—M31 separation di also had to 
be adjusted. In this way, we found that the final mass had 
converged fairly well with just two iterations (Figure 2). 

The changing mass of the MW and M31 meant that one 
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could not trivially convert the separation history d(t) into 
MW and M31 positions (2 /mw(^) respectively). 

However, the instantaneous acceleration of the MW must 
be due to the gravity of Andromeda.^ Thus, the magnitude 

of this acceleration must be a fraction - of the 

total mutual acceleration. This means that 


y MW 


_ ^M31 (^) _ 

^MW H) “ 1 “ ^M31 (^) 

_ 

^MW H) “ 1 “ ^M31 (^) 


(41) 

(42) 


In practice, we solved Equation 32 to determine d{t). 
We found the change in separation over each time timestep 
d{t + At) — d{t) and apportioned this to the MW and M31 
in inverse proportion to their masses at t + | At. 

Our equations are referred to the frame of reference in 
which the origin corresponds to the initial centre of mass 
position (considering only the MW & M31). This makes our 
reference frame inertial. We do not keep track of how the 
centre of mass moves after our simulations start. 

The initial masses of the MW and M31 imply that they 
must have accreted material in some region prior to the start 
of our simulation. Thus, we do not allow test particles to 
start within a certain excluded region. This is defined by 
an equipotential Uexc^ chosen so as to enclose the correct 
total volume (i.e. Vexc pM,i — initial LG mass). 

The density of matter . at the initial time ti includes 
contributions from both baryonic and dark matter. For most 
parameters, the resulting excluded region is a single region 
encompassing both the MW and M31 rather than distinct 
regions around each galaxy. 

The potential resulting from integrating Equation 38 is 


U 


GM 

> - Ln 

/ ^ rp 

j=MW,M31 


V^l + bj^ — 1 


(43) 


We start our test particles on a grid of plane polar co¬ 
ordinates. At some particular angle 0, we consider a sequence 
of trajectories which start further and further out. Trajec¬ 
tories are skipped if they start within the ‘exclusion zone’ 
(U < Uexc at t = ti). Once we obtain a trajectory that fin¬ 
ishes further than 2.15 Mpc from the LG barycentre, we skip 
3 out of every 4 steps as the velocity field is fairly smooth at 
such large distances (Figure 3). Once we reach beyond 3.2 
Mpc, we move on to the next value of 0. This is because we 
do not need the velocity field further than ~ 3 Mpc from the 
LG as there are no target galaxies further away.^ 

We use a fourth-order Runge-Kutta algorithm with an 
adaptive timestep designed to be 30—70 times shorter than 
the instantaneous dynamical time tdyn- This is estimated by 
dividing the distance to each galaxy by the speed of a test 
particle with zero total energy, ignoring the presence of the 
other galaxy. Faced with two estimates of tdyn, we use the 
shorter one in order to maximize the resolution. 

The worst time resolution we use is of the total 



c 


fD . 4 !_I_^^^_I 

^ -1 0 1 2 3 4 

^ Distance from MW-M31 line (x), Mpc 



Figure 3. Top: Local Group velocity field for the case = 0.3, 
Mi = 4 X 10^^Mq and no Centaurus A. Locations of indicated 
galaxies are shown. The MW is just above the origin. Only parti¬ 
cles starting at x ^ 0 (and thus Vx ^ 0) were considered. Thus, the 
presence of particles at a: < 0 indicates intersecting trajectories 
and a disturbed velocity field. Bottom: Radial velocities of test 
particles with respect to the LG barycentre. Vertical lines rep¬ 
resent the distances of M31 and the MW from there. Increased 
velocity dispersion near these galaxies is apparent. Black dots 
on the a:-axis show the distances of target galaxies from the LG 
barycentre. Without proper motions, observations can’t be put on 
such a Hubble diagram as the MW is not at the LG barycentre. 


duration (~13.5 Gyr). This was sufficient for distances 
from each galaxy. At smaller distances, we found that tdyn oc 
r 4 is a good approximation. If required, we improve the time 
resolution in powers of 2 up to a maximum of 5 times (for 
a 2^ = 32x reduction in At). This should provide adequate 
resolution for distances from the MW and M31 greater than 
their respective ‘accretion radii’ Vacc, which we chose to be 
a few disc scale-lengths (Table 2). 


^ This is not strictly true at early times due to the term, but 
we do not expect either galaxy to have accreted much mass at 
that stage because no test particle starts very close to the MW 
or M31. Without mass accretion, the ratio of this term between 
the two galaxies is also inverse to that of their masses. 

^ This requires trajectories starting out to distances of ~0.5 Mpc. 


2.2.2 Including Centaurus A 

Although none of our target galaxies are too close to any 
of the perturbers listed in Table 1 (due to pre-selection by 
Penarrubia et al. 2014), we were still concerned that their 
gravity might have noticeably affected our target galaxies. 
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Name 

b 

1 

dMW 

M 

Gentaurus A 

19.4173° 

309.5159° 

3.8 

4 

M81 

40.9001° 

142.0918° 

3.6 

1.03 

IG 342 

10.5799° 

138.1726° 

3.45 

1.76 


Table 1. Properties of mass concentrations outside the Local 
Group which we considered. Distances are in Mpc and sky posi¬ 
tions are in the Galactic system. The estimate for Gen A is from 
Harris et al. (2010), that of M81 is from Gerke et al. (2011) and 
that of IG 342 is from Wu et al. (2014). Masses are in units of 
10^^ Mq and were obtained from Karachentsev (2005) for Gen 
A and IG 342. For M81, we used Karachentsev & Kashibadze 
(2006). 


To test this scenario, we decided to directly include the most 
massive perturber, Centaurus A. 

Due to the large distance of Cen A from the LG 
(~ 4 Mpc), any peculiar velocity it has is likely to be much 
smaller than its radial velocity. Indeed, this is borne out ob- 
servationally for motion along our line of sight (Karachent¬ 
sev et al. 2007). As a result, Cen A is probably on an almost 
radial orbit with respect to the LG barycentre. Fortunately, 
Cen A is currently located almost directly opposite M31 on 
our sky (cos^ = —0.99, where 0 is the angle on the sky 
between M31 and Cen A). This allowed us to continue using 
our axisymmetric model. 

To initialize each simulation, we need trajectories for the 
MW, M31 and Cen A that match the presently observed dis¬ 
tances to M31 and Cen A. This is done using a 2D Newton- 
Raphson algorithm^ on the initial relative positions of all 
three galaxies along a line. As before, initial velocities were 
found using Equation 37. 

Test particle trajectories were then solved in the usual 
way, with the grid of initial positions centred on the initial 
barycentre of the MW & M31 as before. Including Cen A, 
Equation 40 becomes 


r 


a 

-r — 
a 


E 

j = MW, 
M31, Cen A 


GMj{r-r,) 
r -r.\2 + rs .^)^\r -r.\^ 


(44) 


Eor simplicity, we keep the mass of Cen A fixed at 
4 X 10^^M© (Karachentsev 2005) but still allow the MW 
and M31 to accrete mass. 


2.3 Observations Sample Selection 

Our dataset comes mostly from the catalogue of LG galax¬ 
ies compiled by McConnachie (2012). We used the subset 
of these that were used for a timing argument analysis by 
Peharrubia et al. (2014). This implicitly applies a number 
of criteria. The basic idea behind them was to ensure that 
gravity from the MW and M31 dominates over gravity from 
anything else. Eor this reason, targets > 3 Mpc from the LG 
were not considered. The authors also avoided galaxies too 
close to any major mass concentrations outside the LG. The 
perturbers they considered are listed in Table 1. 


Very close to the MW or M31, there are crossing trajec¬ 
tories and so the model-predicted velocity in such regions is 
not well-defined (top panel of Eigure 3). Eurther away, this 
issue does not arise. Thus, Pefiarrubia et al. (2014) restricted 
their sample to non-satellite galaxies. We further restricted 
their sample by excluding Andromeda XVIII as it is in the 
disturbed region around M31, even if it is unbound. 

We treated HIZSS 3A V B as one object as they are 
almost certainly a binary system. Naturally, we used the 
velocity of its centre of mass, assuming a mass ratio of 13:1 
(Begum et al. 2005). To allow for uncertainty in this ratio, we 
inflated the error on the heliocentric radial velocity (HRV) to 
3.5 km/s. This decision turns out not to matter very much 
because the uncertainty in its distance has a much larger 
effect than uncertainty in its HRV (this is true for most of 
our targets — distances are harder to measure accurately). 

In regions close to the MW and M31, the presence of 
crossing trajectories makes it impossible to uniquely predict 
the velocity of a target galaxy based solely on its position 
(Eigure 3). In such cases, we should reject the target (i.e. 
not use it in our analysis). In practice, we accepted ah of 
our targets in ah cases. We checked the velocity held to 
ensure none of our target galaxies fell in regions with cross¬ 
ing trajectories. Although none of them did so, NGC 3109 
and Antila came close. We tried raising and altering the 
distances to these galaxies within their uncertainties, but 
we still could not get them in a region of crossing trajecto¬ 
ries. In any case, excluding them would not much affect our 
conclusions, as will become apparent later (Eigure 15). As 
a hnal check, IB looked at ah the <7pos and GRV maps (like 
those in Eigure 5) and conhrmed that they were smooth. 

If we had been less fortunate regarding the locations 
of our target galaxies, then we might have rejected some of 
them in some simulations using criteria designed to search 
for intersecting trajectories. The best options seem to be a 
high density of test particles near the present position of 
the target and a high velocity dispersion at that position. In 
this case, we might have to alter Equation 57 by multiplying 
the hrst term on the right by where N is the number 
of target galaxies ‘accepted’ by the algorithm. Additional 
care would have to be taken to ensure the analysis remained 
valid despite N varying with the model parameters (i.e. some 
models might be constrained using fewer observations than 
others).^ 

To convert observations into the same co-ordinates as 
our simulations, we hrst dehned Cartesian xy co-ordinates 
centred on the LG barycentre, with y towards the MW. 
The positions of observed galaxies were converted into this 
system using the equations 

X = (45) 

Vrel = y~yMW ~ ^MW {^MW ' MW^ (^ 6 ) 

is the present distance of the MW from the initial 
position of the LG barycentre. y = i^ direction 

from M31 towards the MW. This is just the opposite of 
the direction in which we observe Andromeda, is the 


^ For stability, we under-relaxed the algorithm, meaning that in 
each iteration, we altered the parameters by 80% of what the 
algorithm would normally have altered them by 


^ If a target galaxy is problematic in only some parts of parameter 
space, then one can simply avoid including it in the analysis alto¬ 
gether, thereby avoiding issues due to N being model-dependent. 
However, this makes poorer use of the available information. 
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distance from the MW to the target galaxy. This is essen¬ 
tially equivalent to its heliocentric distance. We neglected 
the difference that arises because the Sun is not at the cen¬ 
tre of the Milky Way.^ For this reason, we can approximate 
the direction between the MW centre and the target galaxy 
as the direction in which we observe it. 

Although the position of the Sun with respect to the 
Galactic Centre is unimportant for this work, its velocity 
relative to the MW is very important because this velocity 
is ~250 km/s (McMillan 2011). For observational reasons, 
we split this velocity into two components. The MW is a 
disc galaxy, so most of the Sun’s velocity is just ordered 
circular motion within the disc plane. In the absence of non¬ 
circular motions, its speed would be r’c,©- This is known as 
the Local Standard of Rest (LSR) because particles moving 
tangentially at this speed would be at rest in a rotating 
reference frame. 

We temporarily define a 3D Cartesian co-ordinate sys¬ 
tem with X pointing from the Sun towards the Galactic 
Centre, z pointing towards the North Galactic Pole and y 
chosen so as to make the system right-handed. Fortunately, 
y points along the direction of rotation. In this system, the 
velocity of the Sun with respect to the Milky Way (including 
its non-circular motion) is 




Uq 

Vq + Vc,Q 

Wq 


(47) 


The direction towards another galaxy can be deter¬ 
mined from its Galactic co-ordinates using 


d 


MW 


cosb cost 
cosb sin/ 
sin 6 


(48) 


where b is the Galactic latitude and I is the Galactic 
longitude, whose zero point is the direction towards the 
Galactic Centre. Galactic co-ordinates are actually heliocen¬ 
tric, though the distinction is unimportant for very distant 
objects. 

Without proper motions of LG galaxies, only their GRV 
can be constrained. Thus, we project the velocity of the 
Sun with respect to the MW onto the direction towards 
the desired galaxy. This is then added onto its observed 
heliocentric radial velocity. 

GRVobs — HRVobs + Vq • (49) 


This estimate of GRVobs is dependent on the model used 
for Vq, in particular the adopted LSR speed. Thus, for a 
range of plausible values of I’c,©, we stored the resulting 
values of Vq • for each target galaxy. This quantity is 

the difference between its GRV and its HRV. 


2.4 Comparing simulations with observations 

Our simulations yield a velocity field for the LG. To deter¬ 
mine the model-predicted GRV of an observed galaxy, we 
need a test particle landing at exactly the same position. To 
achieve this, we started with whichever test particle landed 


^ Target galaxies are >800 kpc away while the Sun is only ~8kpc 
from the Galactic Centre, well below typical distance errors. 


closest to the targeted final position. We then used a 2D 
Newton-Raphson algorithm on the initial position of this 
particle. The dependence of its final position on its initial one 
was found using finite differencing. For this, we used trajec- 
tories starting at (Xq-J/o). (*o +c**o>2/o) and (aj^.j/o +dj/o). 
with c/Xq = dy^ = 307 pc. Note that we have reverted to the 
usual xy co-ordinates, with y pointing from M31 towards 
the MW and x orthogonal to this direction. 

We considered the Newton-Raphson algorithm to have 
converged once the error in the final position (x, y) was be¬ 
low 0.001% of the distance between the target and the LG 
barycentre. The final velocity of this trajectory {vx^Vy) was 
used to determine the model-predicted GRV of the target 
galaxy. We then corrected this for the motion of the Sun 
with respect to the MW to obtain its model-predicted HRV. 


GRVm 


H RVraodel 


VxX + (vy - y^^){y - 

+ (y - VMwf 

GRVmodel Vq • dj^y^ 


(50) 

(51) 


If the MW or M31 mass is altered, then another sim¬ 
ulation is required. However, if we only wish to alter the 
adopted I’c,©, then this is not necessary. We simply use the 
same GRVmodel but a different Vq (Equation 47). In general, 
this alters HRVmodel- 

To account for distance uncertainties, the target was 
moved to the la upper limit of its observed distance d^yy 
(using the la lower limit instead had a negligible impact 
on our analysis). The Newton-Raphson procedure was then 
repeated targeting the revised position. Once this converged, 
we extracted the GRV from the final trajectory. We took 
the difference between these GRV estimates and called this 
apos- This is the uncertainty in the model-predicted GRV of 
a target galaxy due to uncertainty in its position. 


G'pos — \GRVmodel (^mw GRVmodel (^mw)I 


Here, aa^w uncertainty in the distance to a target 

galaxy. We assume negligible uncertainty in the direction 
towards it, constraining its position to be along a line. The 
velocity field is treated as linear over the part of this line 
where the target galaxy is likely be. Thus, assuming distance 
errors to be Gaussian, GRV model would also have a Gaussian 
distribution. 

To determine apos for M31, we use a slightly different 
procedure because it is not massless. Once we have the time 
history of the MW—Gen A separation, we keep this fixed 
and vary the initial MW—M31 separation to target a revised 
final value. For consistency, we also do not change M(t). The 
effect on the final GRV of M31 is used for apos- 

We expect this procedure to be approximately correct 
because Gen A only affects the GRV of M31 by ~10 km/s, 
making it not crucial to handle tides from Gen A very ac¬ 
curately. It would be possible to do so by recalculating tra¬ 
jectories for all three galaxies with revised target positions, 
but due to numerical difficulties this would probably have 
been less precise. It will become clear later that our results 
are not much affected by the value of apos for M31. 

Altering the MW—M31 separation changes the gravita¬ 
tional field in the rest of the LG, affecting GRVs of objects 
within it. We expect this to be a very small effect and so we 
neglected it. 

We conducted simulations across a wide range of total 
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initial masses Mi and mass fractions in the MW (see 
Table 2). For situations with we took advantage of a 

symmetry that arises between situations with q^ ^ 1 — q^. 
Essentially, the behaviour of M31 in the low-^^ case is equiv¬ 
alent to the behaviour of the MW in the high-g^ case. Thus, 
we did not repeat all our calculations for the latter. 

The positions of observed galaxies were altered in the 
following way: 


X ^ X (unaltered) (53) 

y = Vmsi - Vrei instead of + j/rei) (54) 

where i/rei is still obtained using Equation 46 and is 
therefore unchanged. Equation 54 also applies to Cen A, so 
its final position is now different. This meant we had to find a 
new solution for the trajectories of the MW, M31 and Cen A 
respecting the revised constraint on Cen A. Once this was 
done, we had to deal with altered positions {x, y) for our 
target galaxies by finding new test particle trajectories with 
the right final positions. The final GRVs of these trajectories 
were obtained using Equation 50, but referred to M31 rather 
than to the MW. 

The step we did not repeat was the calculation of the 
LG velocity field. This meant we had a much poorer guess 
for the initial position of each target galaxy. Eor this, we 
simply re-used the values of x and yrei at the initial time in 
the low-^i case. Despite this, our algorithm still converged. 

This procedure implicitly assumes that the accretion 
radii of the MW and M31 are swapped (i.e. that the galaxy 
with the higher mass always has the larger accretion ra¬ 
dius). However, with the very low amounts of mass accreted 
by these galaxies (Eigure 6), this should hardly affect our 
results. This is especially true when considering that our 
analysis tends to disfavour > | (Eigure 7). 

As well as uncertainties due to position (apos) and mea¬ 
surement error on the radial velocity (cr^,^), we also included 
an extra variance term, a extra- This was to account for ef¬ 
fects not handled in our algorithm, for instance interactions 
between LG dwarf galaxies and tides raised by large-scale 
structures (LSS). <7extra is a measure of how much model- 
predicted and actual radial velocities disagree. Including it, 
the contribution to of any particular galaxy i is 

2 (HRV^odel-HRVobsV , 

Xi = I - - - I where (55) 

<7 = \/ 4 “ 4 “ ^extra^ ( 56 ) 


The uncertainty on the motion of the Sun can introduce 
systematic errors into our analysis. Thus, we treated Vc,q as 
another model parameter. However, it is independently con¬ 
strained (239±5 km/s, McMillan 2011). This was accounted 
for using a Gaussian prior, or equivalently by adding an extra 
contribution to Therefore, the total foi* any particular 
model (= combination of model parameters) is 


2 

X = 


E 

Target 

galaxies 


2 I 

Xi + 


^c,0 '^CjO,nominal \ 

^^c,0 / 


(57) 


Models with higher cFextra will necessarily achieve a 
lower . Thus, we can’t use alone to decide which models 
are best. We made use of the fact that the probability of a 



Figure 4. The overall effect on the present GRV of Andromeda 
due to a 10 km/s impulse to its GRV in the past, with the present 
distance to Andromeda constrained. This constraint is maintained 
by altering the initial MW—M31 separation (see text). Because 
of this, impulses applied longer ago have a smaller net effect on 
present motions. 


model matching an individual observation 

1 

P(Observation of galaxy i\ Model) oc — e (58) 


Thus, we recorded both Xi^ ^ind for each observed 
galaxy. The relative model likelihoods were then found using 


P(Model I Observations) oc 



(59) 


If model-predicted and observed HRVs often disagree by 
much more than observational errors, then non-zero values of 
<7extra will bc preferred. Once becomes comparable to the 
number of target galaxies (32), increasing a extra further will 
not much reduce x^- As a result, instead of P increasing with 
(7extra, It will actually start to decrease because of the factors 
of — in Equation 59. One can imagine this as penalizing 

^ i 

models where x^ is so small that such good agreement with 
observations is Too good to be true’. 

In this way, we hoped to constrain a extra- If model- 
predicted and observed HRVs agree well given observational 
uncertainties, then the posterior distribution of a extra would 
peak at or near 0. If that does not occur, then this might 
indicate underestimated observational errors or a failure of 
the model. 

Physically, we expect the main source of astrophysical 
noise contributing to cFextra to be interactions between LG 
dwarf galaxies. However, Andromeda is much heavier than 
them, suggesting that it should be treated somewhat dif¬ 
ferently. This is because a minor merger would affect its 
velocity very little. Thus, whatever the adopted value of 
7 extra for othcr LG galaxies, a smaller value of a extra, M 3 i 
should be adopted for M31. We used =0.1. This 

alters Equation 56 for M31 and thus its contribution to x^. 

We considered the effect of a minor merger with An¬ 
dromeda or the Milky Way in the past. This was modelled 
as an impulse, meaning that we instantaneously altered the 
GRV of M31 at some time in the past. The effect on its 
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present GRV was then determined. For simplicity, Centaurus 
A was omitted and the total LG mass was held constant at 
4 X lO^^M©. This roughly reproduces the present GRV and 
distance of M31. 

One might think that the longer ago the impulse was, 
the bigger its effect on the present GRV of M31, After 
all, pushing the galaxies towards each other increases the 
force between them at later times, further reinforcing the 
original impulse. 

However, this would lead to the constraint on the 
present distance to M31 being violated (in this example, 
it would end up too close). Gonsequently, we had to alter 
the initial separation of the galaxies di compared with a 
non-impulsed trajectory. This tends to counteract the direct 
effect of the impulse. 

The results we obtained for Avj, as a function of the 
impulse time are shown in Figure 4. An impulse applied 
very recently hardly affects df and so di doesn’t have to be 
altered much. Thus, Avj, is almost equal to the impulse. 

For impulses applied longer ago, Avj, rapidly becomes 
very small. The dependence on impulse time is even steeper 
than for Hubble drag (Avj, ^ where a was the cosmic 
scale factor when the impulse was applied). This underlines 
just how difficult it is to alter the present GRV of M31. Gon¬ 
sequently, a realistic model needs to match this constraint 
very well. 

As well as interactions between LG dwarf galaxies, our 
model does not fully account for the presence of large-scale 
structures in the Universe beyond the LG. We attempted to 
include some of these structures in Section 4.3, but others 
remain beyond the scope of this investigation. The lead¬ 
ing order effect of LSS on the LG is to accelerate it as a 
whole without altering the relative velocities between ob¬ 
jects within it. However, LSS also raise tides on the LG, 
affecting the GRVs of our target galaxies. Such effects are 
larger for galaxies further from the MW. As M31 is the clos¬ 
est galaxy in our sample, its GRV should be least affected 
by tides raised by LSS. This further justifies our decision to 
use a value for much smaller than 1. 

extra 


3 RESULTS 

Our analysis works by determining the model-predicted 
GRV of each target galaxy in the LG. A range of models 
are tried out, with different initial LG masses Mi and mass 
fractions in the MW (see Table 2). The results for two 
target galaxies are shown in Figure 5. Each of these GRV 
predictions are converted into a range of HRV predictions 
using Vc,Q within 3a of its most likely value according to the 
work of McMillan (2011). By comparing these with observed 
HRVs, we obtain complementary constraints on the model 
parameters. As we have >> 2 target galaxies, we also test the 
model itself. 

Our simulations allowed the MW and M31 to accrete 
mass. In Figure 6, we show the fraction by which the original 
mass of each galaxy increased. The galaxies only increase 
their mass by a few percent in our simulations. Thus, accre¬ 
tion is unimportant in them. This is mainly because a test 
particle needs to pass within a few disc scale-lengths of the 
MW or M31 for us to consider the particle accreted (Table 
2 ). 


Andromeda grv, km/s 



Sextans A grv, km/s 



0.2 0.4 0.6 0.8 


Figure 5. Simulated galactocentric radial velocities (GRVs) as 
a function of model parameters for two target galaxies. Different 
galaxies constrain a different combination of model parameters, 
with Andromeda mostly telling us the total LG mass and Sextans 
A mostly telling us the fraction of this mass in the MW (because 
the GRV of Sextans A is much higher at low q^). 

This aspect of our models is not totally realistic. If more 
distant approaches were also treated as leading to accretion, 
then the MW V M31 would gain more mass. We do not 
consider this an important effect because we tried a wide 
range of initial masses for both galaxies. 

Figure 7 shows the posterior probability distributions 
of all our model parameters and pairs of parameters based 
on a set of 1128 simulations^ that include Gentaurus A with 
a mass of 4 x lO^^M©. Each simulation was compared with 
observations using 101 values of Vc,q and 201 values of a extra 
(priors are given in Table 2). Of particular importance is the 
posterior on a extras which we constrain to be 45^6 km/s. 
As observational errors are typically ~5—10 km/s and are 
already included in our analysis, this is very surprising. 

We checked if varying the start time of our simulations 
from cii — ^ ^ ^ affected our results. This reduced the 
most likely value of aextra by ~ 1 km/s. Our results are not 
much affected by the epoch at which our simulations are 
started. Some reasons for this are given in Section 4.6). 

^ spanning a linear grid with 24 steps in Mi and 47 steps in q ^, 
although some shortcuts were taken for q^ > \ 
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Figure 6. Fraction of the initial Milky Way mass that it accretes 
by the end of our simulations. Due to symmetry, approximate 
results for Andromeda can be obtained by setting —)■ 1 — . 


We considered a different estimate for the LSR speed 
(238 ±9 km/s, Schonrich 2012). As might be expected, this 
affected a extra by <1 km/s. This is because we consider 
Vc,Q to be well constrained independently of our work. It 
is also apparent that there is very little tension between 
these independent measurements and our timing argument 
analysis (Table 2). 

Our special treatment of M31 forces up a extra to some 
extent as it essentially forces our models to match its GRV 
(given the small uncertainty on r’c,©)- As this may be overly 
restrictive, we redid our analysis using the same value of 
cr extra for M31 as for other LG galaxies (i.e. — 1 

instead of 0.1). This lowered a extra by ~ 2 km/s. 

Our method of handling distance uncertainties is very 
similar to that used by Peharrubia et al. (2014). We rely 
on an assumption that the velocity field is approximately 
linear over the range of positions where the galaxy could 
plausibly be. To test this, we repeated our calculations with 
apos estimated based on how much the simulated GRV of 
each target galaxy changed if we altered its distance from 
its observed value to its Icr lower limit. This gave almost 
identical results to when we used the Icr upper limit instead 
{cr extra dccreascd by ~ 0.1 km/s when using the lower limit). 

Our analysis favours a very low value for q -^, the fraction 
of the LG mass originally in the MW. This is related to 
the fact that observed HRVs tend to systematically exceed 
the predictions of the best-htting model (Figure 9). Thus, 
our analysis will prefer those models that generally lead to 
increased GRVs. Reducing q^ has this effect because it causes 
particles projected orthogonally to the MW—M31 line; to 
curve towards M31 and away from the MW. It also implies 
a faster motion of the MW relative to the LG barycentre and 
a greater distance from there, enhancing projection effects. 

Most of our target galaxies are in fact roughly orthogo¬ 
nal to the MW—M31 line as perceived from the LG barycen¬ 
tre (Figure 3). This might be why q^ seems to have a strong 
impact on GRVs (bottom panel of Figure 5). Thus, one 
might expect our analysis to prefer very low values of q-^, 
which indeed it does. 

Gertain correlations are apparent between some of our 
model parameters. Because we require our models to accu- 


Name 

Meaning and units 

Prior 

Result 

C extra 

Extra velocity dispersion 
along line of sight, km/s 

0 - 100 

45.115;? 

Mi 

Initial MW + M31 mass, 
trillions of solar masses 

2 - 6.6 

4.1T0.3 


Fraction of MW + M31 
mass initially in the MW 

0.04-0.96 

0.14T0.07 

'^C,© 

Circular speed of MW at 
position of Sun, km/s 

239 ±5 

239.5T4.8 

Fixed parameters 

do 

Distance to M31, kpc 

783 ± 25 


Ho 

Hubble constant at the 
present time, km/s/Mpc 

67.3 


Hm,0 

Present matter density in 

3H ^ 

the Universe -F 

0.315 


tti 

Scale factor of Universe 
at start of simulation 

0.1 


'^acc,MW 

'^acc,M31 

Accretion radius of MW 

Accretion radius of M31 

15,337 parsecs 

21,472 parsecs 

Uq 

Vq 

Wq 

See Equation 47 

See Equation 47 

See Equation 47 

14.1 km/s 
14.6 km/s 
6.9 km/s 



Table 2. Priors and Icr confidence levels on model parameters. 
The latter are far from the boundaries imposed by the former, 
showing that our results are not strongly affected by our priors. 
Due to accretion, the present-day masses of the MW and M31 
are ~5% higher than when the simulations start. We use the mea¬ 
surement of do by McConnachie (2012). Cosmological parameters 
are from Planck Collaboration (2015). We obtained Ve,Q from 
McMillan (2011) and the Sun’s non-circular velocity from Francis 
& Anderson (2014). Uncertainty in the latter is much less than in 
the former. We assume r’c,© is within 3cr of its most likely value. 


lately match the observed HRV of M31, our timing argument 
estimate of the LG mass is quite sensitive to anything which 
affects its predicted HRV. M31 is almost directly ‘ahead’ of 
the Sun in its motion around the MW. Thus, for the same 
GRV of M31, increasing Ve,Q decreases its HRV (Equation 
49). To increase its HRV back up to its observed value, its 
GRV would have to be increased, which is only possible in 
a different model where the retarding effect of gravity is 
smaller (i.e. Mi is lower). 

A lot of our target galaxies have HRVs which exceed 
the predictions of the best-fitting model (Figure 8). This 
means that a lower LG mass fits the data slightly better, 
explaining the correlation between Mi and a extra- For the 
same reason, increasing Ve,Q indirectly improves the fit to 
the data, reducing cfextra slightly. 

Some effects are inevitably not considered in our model. 
If they were included, we might achieve a better fit to the 
observations. We consider some of these effects in the next 
section. We pay special attention to tides from objects be¬ 
yond the Local Group (Section 4.3) and the Large Magel¬ 
lanic Gloud (Section 4.4). 
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Figure 7. Posterior probability distributions of the parameters considered in this work (defined in Table 2). In each figure, we marginalized 
over the parameters not shown. For variables plotted against other variables, we show the contours of the probability density corresponding 
to the 1(7 (68.3%) and 2a (95.4%) confidence levels, as well as the most likely pair of values. Rotate figure 90° clockwise for viewing. 
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4 DISCUSSION 

Our analysis reveals an astrophysical noise in velocities of 
LG galaxies that greatly exceeds observational errors. With 
N — 32 targets, the fractional uncertainty in this extra noise 
<y extra should bc ~ This agrees closely with the 

width of the posterior probability distribution of a extra (Fig¬ 
ure 7). 

We considered several factors which could influence our 
analysis. Perhaps most obviously, the LG contains gravity 
from objects other than the MW and M31. For example, 
the non-satellite LG galaxies that we modelled as test par¬ 
ticles in reality exert gravity on each other. This would lead 
to roughly isotropic and random impulses on them. Gon- 
sidering that our analysis is based solely on line of sight 
velocities, we would need to assume typical impulses of 
~ V3 CT extra ~ 80 km/s. 

However, our target galaxies have typical velocity dis¬ 
persions/rotation speeds of < 15 km/s (e.g. Kirby et al. 
2014). For some impact parameters, these galaxies could 
perhaps impulse each other by twice this while avoiding 
a merger. Thus, the high value of a extra inferred by our 
analysis seems difficult to explain as a result of interactions 
amongst the galaxies we considered. 

Additional inaccuracies in our model may arise from the 
effects of large scale structure. Moreover, even distant en¬ 
counters between LG dwarf galaxies can affect their motion. 
The likely magnitude of such effects can be estimated based 
on more detailed cosmological simulations of the AGDM 
paradigm. Gonsidering analogues of the LG in such simula¬ 
tions, it has been found that the dispersion in radial velocity 
with respect to the LG barycentre at fixed distance from 
there should be ~ 30 km/s (Aragon-Galvo et al. 2011). 

Looking at the bottom panel of Figure 3, it is clear that 
our models do not produce such a large velocity dispersion. 
We seem to get cr^ ~ 10 km/s, though this rises slightly 
to ~ 15 km/s once we include Gentaurus A. Thus, even if 
AGDM were correct, it would be reasonable for our analysis 
to infer a extra ~ 25 km/s. 

In this section, we hope to correct some of our model 
deficiencies and make it a more accurate representation of 
AGDM. Table 3 shows some of the effects we consider and 
a rough idea of their contributions to . Gombining every¬ 
thing in quadrature suggests that the objects we consider 
are sufficient to attain a dispersion in the Hubble flow of 
~20 km/s. Thus, a lot of the ‘scatter’ about the Hubble 
flow found by Aragon-Galvo et al. (2011) arises because the 
LG is not spherically symmetric rather than actually being 
a dispersion in velocities at the same position. Nonetheless, 
another ~20 km/s must come from factors we do not con¬ 
sider. This means that values of a extra much greater than 
~20 km/s would be problematic for AGDM. 

Our best-fitting model has = 0.14. Results from this 
model are compared with observations in Figure 8. However, 
we consider it unrealistic for the MW to have only | as much 
mass as M31. Assuming that the virial mass of a halo scales 
as the cube of its velocity dispersion (Evrard et al. 2008) 
and that the ratio of the latter between the MW and M31 
is III = 1.25 (Garignan et al. 2006; Kafle et al. 2012), we 
see that it is unlikely for M31 to have much more than twice 
as much mass as the MW. We believe the best compromise 
between this argument and the low value of preferred by 


Object 

Gontribution 
to (7^ (km/s) 

Gomments 

Section 

MW, M31 & Gen A 

- 15 

10 at 2 Mpc 

2.2 

IG 342 & M81 

- 5 

Table 6 

4.3.2 

The Great Attractor 

- 10 

Equation 62 

4.3.3 


Table 3. Contributions of various sources to cr^, the radial ve¬ 
locity dispersion with respect to the Local Group barycentre at 
fixed distance from there. The Great Attractor leads to a ~40 
km/s range in radial velocities at 3 Mpc (Equation 62). Com¬ 
bining everything in quadrature, we can account for ~20 km/s 
of the ~30 km/s dispersion in the Hubble flow found by Aragon- 
Galvo et al. (2011). This suggests that our model should represent 
AGDM to an accuracy of ^^20 km/s. 



GRV in best fit AGDM model, km/s 

Figure 8. Comparison between predicted and observed GRVs 
based on the most likely model parameters (q^ = 0.14, Mi = 
4.2 X IO^^Mq, Vc,q = 239 km/s). The line of equality is in blue. 


our timing argument analysis (0.14 ±0.07) is found if we set 
q^ =0.2. Thus, when comparing our model predictions with 
observations (Figure 9 onwards), we use the model parame¬ 
ters which best fit the data but with q^ raised to 0.2. This 
raises a extra by ~ 1 km/s and has only a small impact on 
our results, but should make them more realistic. 

Observed GRVs seem to systematically exceed model 
predictions (Figure 8). We used our most plausible model in¬ 
cluding Gen A to subtract model-predicted radial velocities 
from observed ones, yielding AGRV for each target galaxy. 
We then created a histogram of the resulting AGRVs in 
Figure 9, smoothing each data point over its respective un¬ 
certainty. As before, this includes apos and cr^;^. Because it is 
unclear exactly how to convert heliocentric radial velocities 
into Galactocentric ones, we also added crvc,Q ia quadrature 
to all the uncertainties. 

If one assumes that factors outside our model are just 
as likely to raise GRVs of target galaxies as to reduce 
them, then it should be possible to use the population of 
AGRV < 0 galaxies to gain a good idea of how accurately 
our model represents AGDM. The AGRV < 0 galaxies are 
well described by a Gaussian of width 15 km/s (blue area 
in Figure 9). Most of the mismatch is due to Leo A. To ac¬ 
count for tides raised by IG 342 and M81, its radial velocity 
prediction should be reduced by ~ 5 km/s, making it more 
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Figure 9. Histogram showing observed — predicted GRVs (i.e. 
AGRVs) of our target galaxies using our most plausible model 
(q^ = 0.2 instead of 0.14, other parameters as in Figure 8). 
Each data point was convolved with a Gaussian of width a = 
o'pos‘^ + + (Jv^ Q^. We divided our sample into those with 

AGRV < 0 (blue) and those with AGRV > 0 (red). The area 
corresponding to one galaxy is shown as a red square. A Gaus¬ 
sian of width 15 km/s is overplotted as a short-dashed line. This 
matches the AGRV < 0 subsample quite well, especially when 
Leo A is excluded (long-dashed line). 



Figure 10. Positions of target galaxies are shown in the same way 
as the top panel of Figure 3. Distance uncertainties are indicated 
by a thin line. The size of each marker is directly proportional to 
\AGRV\ for the associated galaxy (name beside marker) based 
on the same model as shown in Figure 9, with colour indicating 
sign (red means positive). Uncertainties in AGRV are roughly 
proportional to that in distance. For Leo P, this causes a 30 km/s 
uncertainty, though typical values are much smaller. We expect 
our model to represent AGDM to an accuracy of ~ 25 km/s (see 
text), roughly the same as AGRV of NGG 55. 


consistent with observations (Table 6). In any case, consid¬ 
ering the AGRV < 0 galaxies suggests that inaccuracies in 
our model probably do not exceed 25 km/s, slightly less than 
the Gjj ~ 30 km/s found by Aragon-Calvo et al. (2011) due 
to our careful modelling. Thus, one might expect a Gaus¬ 
sian of around this width to also describe the distribution of 
AGRVs for galaxies with AGRV > 0. 

However, unlike galaxies with AGRV < 0, those with 
AGRV > 0 are not well described by a 15 km/s Gaussian 
(red area in Figure 9). There appears to be a population 
of AGRV > 0 galaxies which might be described by such 
a Gaussian, but in this case we would need perhaps two 
additional populations to fully account for the observations. 
A possible mechanism for generating these populations is 
described in Section 4.6. 

We tried to see if there was any correlation between the 
position of a target galaxy and it’s associated AGRV. This 
is shown in Figure 10, with the size of the marker for each 
galaxy directly proportional to its AGRV. It is assumed that 
the LG is axisymmetric, so positions are shown using the 
same co-ordinate system as our simulations. The uncertain 
distance to each galaxy is indicated by a thin line. 

The objects with the highest AGRVs tend to be fur¬ 
thest from the MW/M31. This might be a sign that tides 
from objects outside the LG are responsible for the discrep¬ 
ancies. As we already included Gentaurus A in our simula¬ 
tions, we might be seeing the effects of other objects. We 
will investigate some possibilities in Section 4.3. In partic¬ 
ular, we will show that IG 342 and M81 are unlikely to be 
responsible for the discrepancies (Section 4.3.2). This is also 
true of the Great Attractor (Section 4.3.3). An explanation 
for this trend is suggested in Section 4.6. 


4.1 Reduced Local Group Mass 

Gomparison with cosmological simulations suggests that the 
timing argument may overestimate the LG mass (Gonzalez 
et al. 2014). Moreover, observed radial velocities tend to 
be systematically more outwards than in our models. These 
considerations suggest that a lower LG mass could help to 
explain the observations. To test this, we removed the effect 
of gravity altogether and used Equation 1 to predict veloc¬ 
ities. As before, we took the barycentre of the LG as the 
centre of expansion and assumed that Vc,q = 239 km/s. 

The MW was assumed to be going towards this point 
at ~90 km/s, which is reasonable given the observed HRV 
of M31 and a plausible mass ratio between the galaxies. In 
theory, the MW should be going away from the LG barycen¬ 
tre in the absence of gravity. However, using the correct 
MW velocity ensures that velocities with respect to the LG 
barycentre are correctly converted into velocities with re¬ 
spect to us. One slightly unusual consequence of this is that 
the model predicts M31 to have a negative GRV. 

GRV predictions obtained in this way are compared 
with observations in Figure 11. Due to the effect of grav¬ 
ity, observed GRVs tend to be less than in a pure Hubble 
flow. Surprisingly, this is not true for some of our target 
galaxies, especially the DDO objects. Other examples of 
this behaviour have been identified recently (Pawlowski V 
McGaugh 2014a). Thus, reducing Mi does not explain the 
observations, at least if considered on its own. 

Moreover, there is limited scope to alter Mi because it 
is tightly correlated with the present GRV of M31 (Figure 
5). A model needs to match its GRV fairly well because it is 
unlikely that a minor merger with M31 or the MW could 
have substantially affected their relative motion. Even if 
such an event did occur, its net effect on the present GRV of 
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Figure 11. Comparison between observed GRVs and the predic¬ 
tions of Equation 1, with the centre of expansion at the barycentre 
of the LG assuming = 0.2. The MW is taken to be moving 
towards this point at 90 km/s. Notice that some galaxies are 
moving outwards even faster than a pure Hubble flow (blue line). 


6 



Figure 12. This shows a simple ID model for the fractional over¬ 
density S at the start of our simulations, plotted in co-moving co¬ 
ordinates X- The MW and M31 are treated as one object which 
formed from the over-density at low x- The over-density must be 
surrounded by an under-dense region so that the two cancel out 
on large scales (i.e. f S (x) X^^X = 0)- To avoid negative density, 
S- > —1. We also show some test particle trajectories based on 
solving Equation 60. Particles starting further out than B simply 
trace the cosmic expansion (x is constant) because the mass en¬ 
closed interior to their radius is the same as in a homogeneous 
universe. This is not true for particle C — the region it encloses 
is over-dense, so x of particle C decreases slightly. The effect is 
larger for particles starting closer to the LG barycentre (e.g. D). 


M31 would be greatly diminished unless it occurred recently 
(Figure 4). 

In our models, the mass in the LG is present not only 
in the MW and M31. We assume that the rest of the LG 
(RLG) contains a uniform distribution of matter with the 
same density as the mean cosmic density of matter p. At 
this density, a sphere of co-moving radius 2.9 Mpc would 
have a mass equal to that of the MW and M31, assuming 
^MW + M31 = 4 X 10 Mq. 


At early times, it is possible that a small region encom¬ 
passed the material that would later end up in the MW and 
M31. A surrounding under-dense region would be required 
so that the mean density in the union of both regions was 
p. This is depicted schematically in Figure 12. 

If we assume the under-dense region was completely 
empty (i.e. S- = —1), then it would have to extend out 
to a co-moving radius of 2.9 Mpc. As a result, there would 
be no mass in the RLG, assuming this was defined to have a 
radius below 2.9 Mpc. It can be seen from Figure 3 that all 
our target galaxies have distances from the LG barycentre 
of < 2.9 Mpc. Thus, they could all be in a void. 

However, one must bear in mind that test particles are 
retarded by the gravity of the MW and M31. This reduces 
the co-moving volume spanned by a cloud of test particles. 
Thus, if the RLG is not completely empty, then the LG con¬ 
tains material initially outside its present co-moving volume. 

To investigate the interplay between these effects, we 
now consider the opposite limit in which |(5_| ^ 1. This 
corresponds to a much larger under-dense region surround¬ 
ing the MW and M31 at the start of the simulations (~500 
million years after the Big Bang). To better understand this 
case, we solved some test particle trajectories assuming a 
point mass in an otherwise homogeneous universe. We kept 
fixed the mass enclosed interior to the radius of any given 
test particle. This makes the equation of motion^ 

^ = -^^+Ho"^A.or (60) 

4:71 Q 

Meff = M-h ^r.^p. (1-h S-) (note: (5- < 0) (61) 

O 

p. is the mean density of matter in the Universe at the 
time our simulations are started. Me// includes both the 
point mass M and any material originally present at radii 
below the initial radius of the test particle. We assume that 
Meff remains constant because the system avoids crossing 
of shells. This can be achieved if the massive object accretes 
any objects that come sufficiently close to it, rather than 
just letting them escape on the other side. 

To obtain a final distance from the LG of 2.9 Mpc, we 
need an initial distance of 0.46 Mpc for a starting time cor¬ 
responding to when a. =0.1. This means that 16 x 10^^Mq 
would end up within the RLG at the present time. If the 
RLG were to contain matter at p, then it would only contain 
4 X 10^^M©. Thus, the RLG might have up to = 4 

times as much material as was assumed in our calculations. 

It is difficult to know how much mass is actually present 
in the RLG. There might be a diffuse component of dark 
matter or concentrations of it that have no detectable stars. 
Some regions are difficult to survey because of e.g. the disc 
of our Galaxy. Recently, significant amounts of hot gas have 
been discovered around the MW (Salem et al. 2015) and 
around M31 (Lehner et al. 2015). 

For these reasons, we assumed neither of the extreme 
cases just outlined. Instead, we used an intermediate as¬ 
sumption that the RLG contains matter with a mean den¬ 
sity of p and little density variation. Roughly speaking, this 
corresponds to S- = —0.43 and an under-density out to 3.5 


^ It can be verified that a pure Hubble flow is recovered for the 
case M = S- = 0. 


MNRAS 000, 1-27 (2016) 









Dynamical History of the Local Group in ACDM 17 


Mpc. We think this is reasonable considering the distances 
to major mass concentrations just outside the LG (Table 1). 

We investigated whether altering this assumption might 
affect our conclusions regarding the inferred value of cr extra- 
To do this, we assumed the extreme case that the RLG 
has no mass. This means that the ^ term present in the 
equations of motion (e.g. Equation 32) should be replaced 
with We re-ran our entire analysis using equations 

of motion altered in this way. 

If we used the same procedure as before to prevent test 
particles starting too close to the MW/M31, then we would 
end up with no test particles within ~2.9 Mpc of the LG 
barycentre. In this case, there would be no way to obtain 
HRV predictions for our target galaxies, consistent with the 
assumption of an empty RLG. Glearly, this assumption is 
wrong at some level. Thus, we allow test particle trajectories 
to start anywhere as long as they end up at the correct 
position. 

Altering the setup of our simulations in this way re¬ 
duced a extra by ~7 km/s. Individual radial velocities are 
often increased by larger amounts. This tends to reduce the 
discrepancy with observations. However, the overall effect is 
small because, for the same MW and M31 mass, the GRV 
of M31 is increased. To bring it back down to the observed 
value, the mass of the MW and M31 have to be increased, 
reducing the predicted GRVs of other LG galaxies. 

We checked this explicitly by comparing the 
marginalised posterior probability distribution of Mi. 
Assuming the RLG has a mean density of p, the total initial 
LG mass in units of 10^^Mq is 4.2 ± 0.4. However, if we 
assume an empty RLG, this rises to 5.2 ± 0.4. This seems 
rather high, but is in line with similar calculations by other 
workers (Partridge et al. 2013).^ We suggest that this result 
points towards a RLG that can’t be considered empty 
for the purposes of the timing argument. However, more 
reasonable values for Mi are obtained if one includes the 
kinematic effect of a sufficiently massive Large Magellanic 
Gloud (Figure 18). 

4.2 Increased Hubble constant 

Another way to increase model-predicted HRVs is to in¬ 
crease Hq. The cosmological value seems to be fairly well 
constrained (Planck Gollaboration 2015). Once certain bi¬ 
ases are taken into account, this measurement seems to be 
consistent with surveys of Type la supernovae (Rigault et al. 
2015). However, there is also some cosmic variance: under- 
dense regions of the Universe expand faster than the average. 
If we are in such a region, this would lead to the value of Hq 
appropriate for the local Universe being higher than that for 
the Universe as a whole (Wojtak et al. 2014). 

To account for this possibility, we performed another 
simulation with raised by 5 km/s/Mpc. However, we were 
careful to bear in mind that Planck gives a tight constraint 
on the age of the Universe. To avoid altering this, we had 
to further adjust the adopted cosmology. For simplicity, we 
kept this flat. The parameters used are shown in Table 4. 

The resulting posterior on a extra is shown in Figure 

^ Part of the difference arises because Cen A is not usually in¬ 
cluded in timing argument analyses of the LG. 


Parameter 

Old value 

New value 

Ho 

67.3 km/s/Mpc 

72.3 km/s/Mpc 


0.315 

0.243 

Ua,o 

0.685 

0.757 


Table 4. Alterations to cosmological parameters in Table 2 for 
the high Hq model shown in Figure 13. In both cases, the universe 
is flat and equally old (13.81 Gyr). 



Figure 13. The posterior on a extra is shown for two plausible 
values of the Hubble constant Hq. The age of the universe is 
the same in both models, with parameters adjusted accordingly 
(Table 4). Raising Hq by 5 km/s/Mpc reduces cr extra by ~5 km/s. 

13. As might be expected, increasing Hq lowers the inferred 
value of (7extra, but Only by ~5 km/s. This is similar to the 
effect of assuming the rest of the LG is empty instead of 
filled with matter at a density of p (Section 4.1). This is 
reassuring as the simulations work in slightly different ways. 

Of course, it is only possible to count this reduction in 
(7extra oncci to account for the RLG being less dense than in 
our models, one can either alter the equations of motion to 
make the RLG empty or one can raise Hq slightly. Whichever 
method one prefers to use, the effect is not sufficient to ex¬ 
plain the observations, although it does help. 

One thing that may be in favour of models with an 
under-dense RLG is the inferred value of , the fraction of 
the LG mass in the MW. When we tried to make the RLG 
empty by altering the equations of motion (Section 4.1), 
our analysis preferred 0.381 q q 5 instead of 0.14 ± 0.07. One 
might expect a similar effect to occur when we raise Hq — 
after all, the effect on cFextra is very similar. However, our 
calculations show almost no change in the inferred value of 
due to a higher Hubble constant. 

4.3 Tides from objects outside the Local Group 

4-3.1 Centaurus A 

To better understand the effect of Gen A on our results, we 
repeated our analysis without including it. The results are 
shown in Table 5. Broadly speaking, the results are similar 
in both cases, although there are some subtle differences. 

Being close to the MW—M31 fine, Gen A pulls the MW 
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Parameter Prior Posterior Posterior with 

& units without Cen A Centaurus A 


oextra, km/s 

0 - 100 

54+8.9 

46+'^-4 

Mi, IQi^M© 

2 - 6.6 

q ^9+0.35 

4 14+0-35 


0.2 - 0.8 

0 20+0 060 

0 

to 

0 

1 + 

0 0 

b 

CD 

Uc,o, km/s 

239 ±5 

242.0I4;® 

239.91^;® 


Table 5. Effect of Centaurus A on posteriors of model parame¬ 
ters. Both analyses shown here have a uniform prior on q-^ over 
the range 0.2—0.8. 



Figure 14. Posterior probability distribution of Cextra under dif¬ 
ferent model assumptions. Including Centaurus A reduces a extra 
by ~8 km/s. Also raising Hq as described in Table 4 reduces 
O'extra by a further km/s. 

and M31 apart, increasing the GRV of M31 by ~10 km/s. 
To bring it back down to the observed value, Mi would 
need to be increased by ~ (Figure 5). This is in¬ 

deed roughly what happens to the posterior on Mi. Other 
effects are harder to understand, such as why including Cen 
A leads to better agreement with the LSR speed measured 
by McMillan (2011). 

The posterior distribution of a extra is shown in Fig¬ 
ure 14. Including Cen A reduces its most likely value from 
54^7 km/s to 46^6 km/s. If Hq is also increased to 72.3 
km/s/Mpc, then Oextra is further reduced to dll® km/s. 

We only tried one possible mass of Cen A (4 x 10^^M©). 
Including it at this mass reduces a extra by ~ 8 km/s. It is 
possible that adopting a higher mass would reduce it fur¬ 
ther.^ However, it is unlikely that Cen A is more massive 
than 5x lO^^M© (see Figure 1 of Karachentsev (2005)). This 
is only 25% higher than our adopted mass. Thus, using the 
highest plausible Cen A mass rather than our adopted value 
might well reduce a extra, but only by another ~ 2 km/s. 

The inclusion of Centaurus A affects the Hubble dia¬ 
gram for the LC, increasing This is a tidal effect and 
is therefore larger at greater distances. At 3 Mpc, we found 
a range in radial velocity from the LC barycentre of ~ 70 

^ Though it might not, see Figure 17. 


km/s, falling to perhaps half that at 2 Mpc. This corre¬ 
sponds to cr^ ~ 10 — 15 km/s. Although this is below the 
~30 km/s found in cosmological simulations (Aragon-Calvo 
et al. 2011), it does suggest that some of the ‘scatter’ about 
the Hubble flow found in such simulations can be accounted 
for using an axisymmetric model rather than a spherically 
symmetric one. 


4.3.2 IC342 andM81 

Including Centaurus A improves the fit to the data slightly 
but still leaves a very poor ht. As it is the most massive 
perturber, this suggests that tides can’t explain the dis¬ 
crepant HRVs. To check this conclusion, we cross-correlated 
the discrepancies in the HRVs with the distances between 
our target galaxies and the remaining perturbers in Table 
1. This is shown in Figure 15. The discrepancy seems to be 
larger for objects closer to IC 342 or to M81. Thus, we tried 
to see if tides from these objects might help to explain the 
observations. 

We provide two ways of estimating the effects of tides 
raised by IC 342 and M81 on the Local Croup. First, we 
treat each perturber as the only object in the universe. We 
solve test particle trajectories in the usual way and target 
a particular hnal separation with the perturber. We then 
record the peculiar velocity of this trajectory. Using the 
perturber masses in Table 1, the results obtained in this 
way are indicated in km/s on the gridlines of Figure 15. 

In this very simplistic model, the effect of each perturber 
is just an extra velocity towards it with the calculated mag¬ 
nitude. However, the direction of this velocity is not directly 
away from the MW. For example, IC 342 should hardly affect 
the CRV of KKH 98 because, as perceived by KKH 98, the 
MW and IC 342 are almost at right angles (angle ~ 89.4°). 

The perturber would also have a small effect on the mo¬ 
tion of the MW, this being ~15 km/s towards each perturber 
in the context of this model. For a target near a perturber, 
one expects them to also be nearby on the sky. Thus, the 
MW would be pulled towards the target to some extent, 
reducing its CRV. This might be why our more detailed 
model for tides (see below) often predicts that they would 
reduce the CRVs of target galaxies. 

Our more detailed model involves two gravitating 
masses. We treat the MW V M31 as a single object with 
mass 4 x 10^^M© and assume — 0.2. This object repre¬ 
sents the LC. We put the LC and the perturber along the y- 
axis and solve both objects forwards using Equation 32 (the 
relevant mass is that of the MW+M31+perturber). Their 
initial separation is varied so as to get a hnal separation 
equal to the observed distance between the LC barycentre 
and the perturber. 

We then determine how a target galaxy would ht into 
this picture. We solve a test particle trajectory so that it 
ends up at the correct distance from the LC barycentre and 
at the correct angle to the perturber as perceived at the LC 
particle.^ The hnal CRV of the test particle is determined 
using Equation 50, referred to the LC particle rather than 
the MW. 


^ A 2D model is sufficient for this as there are three particles. 
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Figure 15. Observed — predicted HRVs (i.e. AHRVs) of indi¬ 
cated galaxies as a function of distance from the perturbers in 
Table 1. The radius of each marker oc the discrepancy, which 
we list below the name of each galaxy. Blue indicates a mea¬ 
sured HRV below that in the most plausible model = 0.2, 
Mi = 4.2 X 10^^ Mq), while red shows the opposite. The num¬ 
bers on the gridlines show the peculiar velocity in km/s at that 
distance from each perturber if it was the only object in the 
universe (see text). A more detailed model for how perturbers 
affect observations is shown in Table 6. 


Galaxy 

AHRV 

(km/s) 

Effect 
of M81 

Effect 
of IC 342 

DDO 99 

HOT 13 

-6.5 

-10.3 

DDO 190 

104 ±5 

-7.5 

-9.8 

KKH 98 

100 ± 12 

-3.0 

-9.5 

DDO 125 

95.5 ±4.5 

-5.9 

-10.6 

NGC 404 

89 ±37 

-3.6 

-13.2 

Tucana 

48 ±6 

2.0 

5.4 

NGC 3109 

38.1 ±5.1 

-1.9 

1.4 

ESO 294-GOlO 

-16.8 ±5.8 

3.6 

4.9 

IC 4662 

-24 ± 14 

3.1 

9.4 

IC 5152 

-32.3 ±4.0 

3.7 

7.1 

Leo A 

-46.8 ±4.4 

-2.3 

-3.1 


Table 6. AHRVs for the LG galaxies most discrepant with 
our model. Error budgets are found by adding apos and in 
quadrature. We also give an estimate in km/s for how much M81 
and IC 342 might have affected the GRV of each galaxy (method 
described in text). 


To determine the effect of the perturber, we then (ef¬ 
fectively) reduce the perturber mass to 0 and repeat the 
calculation. The final GRV of the test particle is compared 
between the two simulations. Some results from this proce¬ 
dure are shown in Table 6. 

The combination of large distances from the perturbers 
(> 2 Mpc) and projection effects reduce how much tides 
might have affected the GRVs of target galaxies. As a re¬ 
sult, tides from IG 342 and M81 can’t explain the very high 
HRVs of targets such as the DDO objects in the context 
of this model. In fact, for several galaxies like these, tides 
seem to reduce GRVs and thus make the discrepancy even 
worse. Thus, we do not believe that tides are responsible for 


the discrepancies, assuming we have reasonable perturber 
masses (Table 1) and a good method of estimating their 
effects. 


4-3.3 The Great Attraetor 

There are additional structures in the Universe on a larger 
scale which might be pertinent to our analysis. In particular, 
the Local Group as a whole has a velocity of ~630 km/s 
with respect to the surface of last scattering (Kogut et al. 
1993). It is thought that this is mostly due to the gravity of 
the Great Attractor (GA, Mieske et al. 2005). Assuming a 
distance of 84 Mpc, it is clear that tides raised by the GA 
can have a non-negligible impact on motions within the LG. 

As the GA is much more distant from the MW than ob¬ 
jects within the LG, we use the distant tide approximation. 
Treating the MW and other target galaxies as freely falling 
in the gravitational field of a distant point mass, the change 
in the Galactocentric radial velocity of a target galaxy is 

AGRVga = (3 cos^ e-1) for d < do a (62) 

dcA 

Here, do a is the distance to the GA while 0 is the angle 
on our sky between it and the target galaxy, which is at a 
heliocentric distance d. The GA is assumed to have caused 
the LG to gain a peculiar velocity of Vpec,LG — 630 km/s. 
We take the GA to be in the direction I = 325°, b = —7° in 
Galactic co-ordinates (Kraan-Korteweg 2000). 

For 0 close to 0 or 180°, the GA tends to increase 
GRVs. However, for 0 close to 90°, the GA reduces GRVs. 
This arises because both the MW and the target galaxy fall 
towards the perturber at similar rates. As their co-moving 
distance from the GA decreases, so also does their co-moving 
distance from each other. 

Due to the GA, a test particle started with the same 
initial conditions will end up with an altered position as 
well as velocity. Thus, the initial position must be altered to 
match a hxed hnal position. This reduces the effect of the 
GA on predicted velocities of target galaxies (a similar effect 
is shown in Figure 4). 

Because of Hubble drag, present peculiar velocities are 
mostly sensitive to tides at late times. Thus, we expect Equa¬ 
tion 62 to provide a reasonable approximation as long as we 
have accurate distances to the relevant objects and know 
their sky positions. 

Although we may overestimate the magnitude of 
AGRVga, it is much harder to get its sign wrong. This is 
because the sign is dependent on the factor of (3 cos^ ^ — l), 
a quantity sensitive only to the (usually well-known) sky 
positions of relevant objects but not to their distances. The 
trajectory of a distant LG galaxy is unlikely to have deviated 
much from a radial orbit. The GA is much more distant so it 
was probably always in much the same direction on the sky. 
As a result, we expect 0 to never have been much different 
to its present value. Although this may not be true at very 
early times, Hubble drag makes the system ‘forget’ about 
the forces acting at such times. 

In Figure 16, we show how much tides from the GA 
would likely affect the HRV of each target galaxy. It is in¬ 
teresting to see the results for the galaxies with the highest 
AHRVs, in particular the DDO objects (Table 6). Because 
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Figure 16. AHRVs are plotted against our estimate for how 
much the Great Attractor might have increased the HRV of each 
galaxy. The distance d used here is heliocentric. Notice that the 
DDO objects would likely have their radial velocities reduced by 
tides from the Great Attractor. 


of their positions on the sky, tides from the GA would ac¬ 
tually reduce their GRVs. This makes it more difficult to 
explain their high observed HRVs. 

Perhaps more important is the lack of any apparent 
correlation between HRV discrepancies and the effect of tides 
raised by the GA. This suggests that it can’t reconcile the 
differences between our best-fitting simulation and observa¬ 
tions. In fact, it would probably make matters worse as it 
reduces HRVs for 3 out of the 4 objects which likely have 
AHRV > 90 km/s (i.e. AHRV > 3cr^). 

To estimate how much the GA might affect a extra, we 
adjusted model-predicted GRVs of all our target galaxies 
using Equation 62. We then re-ran our statistical analysis. 
This raised a extra by km/s (or 4 km/s if the RLG is 
assumed empty). The actual effect of the GA is probably 
smaller because it affects final positions of test particles as 
well as their velocities. Still, it seems likely that the GA 
makes it harder rather than easier to explain observed HRVs. 

Although the only large-scale structure we consider ex¬ 
plicitly is the Great Attractor, it is of course possible to per¬ 
form a similar analysis for an external perturber in another 
direction. Indeed, all possible directions can be investigated 
using a grid method. This was done by Peharrubia et al. 
(2016), who used a ID model for the LG (see Section 4.7). 
Fortunately, tides raised by LSS are more important towards 
the edges of the LG. Here, the greater distance from the 
MW and M31 makes it more realistic to consider them as 
a single point mass. Thus, one might expect their analysis 
to be reasonably sensitive to tides raised by LSS. As a re¬ 
sult, it is important to note that they ‘found no statistically 
meaningful deviation between the velocities predicted by the 
point-mass model and the location of galaxies on the sky.’ 


4.4 Kinematic corrections due to massive 
satellites 

Massive satellite galaxies of the MW can affect our timing 
argument analysis because of an indirect kinematic effect. 


Instead of dealing with just the MW, we should really deal 
with the MW system (= MW + satellite). The brightest 
satellite galaxy of the MW is the Large Magellanic Gloud 
(LMG). This is ~ 50 kpc from the Sun (Pietrzyhski et al. 
2013). Being much fainter than the MW, we expect it to be 
much less massive. As a result, it shifts the barycentre of the 
MW system by < 10 kpc. 

Previously, we neglected errors that arise due to the 
heliocentric directions towards other LG galaxies not being 
the same as the directions from the centre of the MW. This is 
because the Sun is only ~ 8 kpc from there (McMillan 2011). 
Similarly, we also neglect any errors that arise due to the 
LMG altering the position of the MW system’s barycentre. 
This is because even the nearest target galaxy is ~ 800 kpc 
away. Moreover, the directions from the Sun towards the 
Galactic Gentre and towards the LMG are almost orthog¬ 
onal, meaning that the errors due to these approximations 
would add in quadrature rather than linearly. 

Unlike the position of the MW system’s barycentre, its 
velocity may be significantly altered by the LMG. Gonse- 
quently, we determined its space velocity with respect 

to the MW. This requires knowledge of its heliocentric radial 
velocity (McGonnachie 2012) and its proper motion (Kalli- 
vayalil et al. 2013) multiplied by its distance (Pietrzyhski 
et al. 2013). The velocity of the Sun with respect to the 
MW is also required (Table 2).^ Using these references, we 
found that the speed of the LMG with respect to the MW 
is 319,845.6 m/s of which 229,401.5 m/s is directed towards 
the North Galactic Pole. Importantly, the component of this 
velocity directly away from M31 is 241,223.4 m/s. 

To apply a kinematic correction for the motion of the 
LMG, we note that the velocity of the Sun with respect to 
the MW should now be altered to its velocity with respect 
to the barycentre of the MW system. 


v© 

Qlmc 


'^O ~ frecoil Qlmc'^LMC 

Mlmc 

Mmw + Mlmc 


(63) 

(64) 


Here, is the velocity of the LMG with respect 

to the MW disc. Note that the MW mass does not 

include a contribution from the LMG mass The — 

sign in Equation 63 arises because we are correcting for the 
recoil induced by the LMG on the MW. If the LMG were 
bound to the MW and the two were orbiting their common 
centre of mass, then a simple application of Newton’s third 
law would show that we should set frecoU = 1- 

ft is possible that the LMG is not bound to the MW but 
rather is on a first infall trajectory (e.g. Besla et al. 2007). 
Indeed, its high speed relative to the MW means it is unlikely 
to be a gravitationally bound satellite galaxy (Wu et al. 
2008). Although the magnitude of is now believed to 

be smaller than the ~380 km/s assumed in this work, other 
considerations continue to suggest that it is unbound (see 
section 6.4 of Kallivayalil et al. (2013)). Thus, most of its 
velocity might have been present even when it was far from 
the MW. In this case, only part of its velocity would have 
been gained due to gravity from the MW. As a result, the 


^ The important quantity here is actually the velocity of the LMG 
with respect to the Sun, so it is essential to use the same vq as 
in the rest of our analysis. 
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recoil of the MW induced by the LMC would be less than if 
the MW and LMC were bound. 

To account for this possibility, we introduce the param¬ 
eter frecoil , the fractiou of the momentum of the LMC that 
has been gained due to gravity from the MW. We assume 
that the direction of the recoil induced by the LMC on the 
MW is aligned with the MW—LMC relative velocity. The 
effect of the LMC on our analysis is maximized if we set 
frecoii = 1 and assume the LMC is bound to the MW. 
The validity of this assumption remains an open question. 
However, to put an upper limit on the effect of the LMC, 
we will make this assumption. 

In applying a kinematic correction for the LMC, another 
major uncertainty is its mass relative to the MW. Recent 
rotation curve measurements of the LMC based on both 
radial velocities and proper motions indicate a ffatline level 
in its outer regions of ~ 90 km/s (Kallivayalil et al. 2013). 
Extrapolating to a tidal radius of 25 kpc as suggested by 
this work, we obtain an enclosed mass of 4.7 x 1O^°M0. The 
actual value is likely to be smaller because other studies 
indicate a slower-rotating LMC (Alves & Nelson 2000). 

The LMC mass can also be estimated using an abun¬ 
dance matching technique. This yields a pre-infall mass of 
~ 1.9 X 10^^ Mq (Boylan-Kolchin et al. 2010). Not all of 
this mass can get as close as 50 kpc to the MW and exert 
a force on it. This is because the outer parts of the LMC’s 
dark matter halo have likely been tidally stripped due to 
its close approach of the MW. The work of Boylan-Kolchin 
et al. (2010) suggests that we should reduce the pre-infall 
mass of the LMC by a factor of ~ 3.5 to account for this 
(see their Figure 11). This makes both estimates of the LMC 
mass agree. 

However, several recent investigations suggest a much 
higher LMC mass, which may be possible if it has not been 
tidally stripped to a signihcant extent. This is tied to the 
issue of whether the LMC is on its first infall into the MW, 
as first suggested by Besla et al. (2007). Those authors 
conducted further investigations into this possibility (Besla 
et al. 2010, 2012). Recently, it was shown that a first infall 
of a massive LMC could induce a recoil on the MW of as 
much as ~70 km/s (Gomez et al. 2015), corresponding to 
Qlmc ^ ^ high LMC mass is also hinted at by the 

discovery of stellar streams around the Magellanic Clouds 
(Belokurov & Koposov 2016) and by its high star formation 
rate, suggestive of a first infall (Tollerud et al. 2011). 

A very massive LMC would exert strong tides on the 
disc of the MW, perhaps warping it more than is observed. 
Assuming a bound LMC, it proved possible to reproduce 
important properties of the observed warp with a fairly low 
LMC mass of just 20 x lO^M© (Weinberg & Blitz 2006). If 
the LMC was instead on its first infall, it would only recently 
have had a substantial effect on the MW. This might be 
compensated by a higher LMC mass. The interplay between 
these effects deserves further investigation. 

To incorporate the LMC into our analysis, we assumed 
that the relevant MW mass for the purposes of the timing 
argument is the combined mass of the MW and the LMC. 
Even if the LMC was quite far from the MW in the past, 
it seems likely that other LG galaxies were much further 
still, so that the MW and LMC can be treated as a single 
point mass. Neglecting the small increase in MW mass due 



Figure 17. The most likely value of a extra and its la uncertainty 
are shown here as a function of the LMC mass (assuming frecoii = 
1, see text). The LMC is included via Equation 63. Notice that 
(yextra eventually increases with the LMC mass because very high 
LMC masses imply a very large kinematic correction for it. 


to accretion (Figure 6), this means that 


Qlmc 


Mlmc 

q^Mi 


(65) 


In models with a very low total LG mass Mi and a 
very small fraction of this in the MW, it is possible to 
get q^Mc > 1- avoid this occurring, we calculated qj^^c 
using Equation 65 and then capped its value at 0.3. This is a 
very generous upper limit on the ratio between the LMC and 
MW masses — high-resolution ACDM simulations indicate 
that it is very unlikely to find a sub-halo with > 10% as 
much mass as the main halo (Boylan-Kolchin et al. 2010). 
Moreover, virial masses scale approximately as the cube of 
rotation velocities (Evrard et al. 2008). Assuming the MW 
rotation curve in the DM-dominated regions is above 180 
km/s (Kaffe et al. 2012), this suggests a LMC mass of below 
~ I that of the MW. 

To check if our imposed upper limit on q^Mc 
fecting our analysis, we determined the best-fitting values of 
Mi and q^ for each value of that we tried. We then 

verified that the most likely total MW system mass (M^^^) 
always comfortably exceeded ^ times the LMC mass (in 
fact, it was never below 4.5M^^^ — see Figure 18). This 
indicates that our analysis should not have been greatly 
affected by our decision to cap q^Mc 

We investigated a wide range of LMC masses 
(0—250x lO^M©) and tried 101 linearly spaced values in this 
range. Each time, we recalculated the probability density 
function over the other 4 model parameters, meaning that 
we essentially added a hfth parameter using a grid method. 
As we are only including the kinematic effect of the LMC, it 
is not necessary to repeat our dynamical simulations. Only 
the statistical analysis had to be redone. The posterior prob¬ 
ability distribution was then marginalized over each of the 
model parameters one at a time to look for trends with the 
LMC mass. 

Probably the most important result of this investiga¬ 
tion is that the overall fit to the observations is not much 
improved. We quantify this using a extras which is reduced 
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Mmw+Mm31 


Figure 18. The red line shows the locus of the most likely values 
of Mi and as a function of LMC mass. Crosses show uncer¬ 
tainties on each parameter for 6 different LMC masses. Figure 7 
suggests that uncertainties in Mi and q-^ are nearly uncorrelated. 


slightly once the LMC is included (Figure 17). However, for 
very large LMC masses, the kinematic correction it induces 
becomes very large, thereby worsening the ht to the data. 
Thus, including the LMC can’t reduce aextra by even as 
much as its formal uncertainty. 

The correction for the LMC is implemented by altering 
Vq according to Equation 63. At given adds a 

constant vector to the predicted velocities of all LG galaxies 
with respect to the Sun. If we set qj^Mc — ^-2, then the 
magnitude of this vector is 64 km/s, comfortably exceed¬ 
ing a extra- Nonetheless, including the LMC hardly reduces 
cr extra- This is because the galaxies we identihed as having 
anomalously high radial velocities (Table 6) are in several 
quite different sky directions. Indeed, we conhrmed that 
assuming a large LMC mass causes some galaxies to have 
HRVs very substantially below the predictions of the best- 
htting model. 

The motion of the LMC with respect to the MW disc 
is mostly along the MW—M31 line, so one expects a strong 
effect on the implied total LG mass Mi. This is clearly borne 
out by our analysis (Figure 18). The lower LG mass resulting 
from including the LMC is now more consistent with the 
works of Peharrubia et al. (2014) and Diaz et al. (2014) — 
both give values around 2.5 x 

Our analysis clearly prefers a non-zero value for qj^^^ 
(Figure 17). Thus, at low LMC masses, the analysis prefers 
low MW masses to force up q^^c towards its preferred 
value. The opposite occurs at high LMC masses — a rapid 
increase in q^ is required to raise the MW mass and hold 
down the kinematic correction due to the LMC. This is be¬ 
cause the ht to the data is worsened if this correction is too 
large. 

The LMC also affects the tangential velocity of M31 
with respect to the MW system. Considering the M31 proper 
motion measurement of van der Marel et al. (2012a), it is 
likely that including the LMC slightly increases the tangen¬ 
tial velocity of M31, though its radial velocity is increased 
far more. Still, M31 should remain on a nearly radial or¬ 
bit with respect to the MW system if the LMC is given a 


frecoil 


0 

239.53 ± 4.82 km/s 

125x109Mo 

239.44 ± 4.72 km/s 

250x109Mo 

238.55l4g7 1^™/® 


Table 7. The most likely value of the LSR speed Vc,q with its 
1(7 uncertainty is given as a function of the assumed LMC mass. 
The prior constraint is 239 ± 5 km/s (McMillan 2011). There is 
no tension between this value and that suggested by our analysis. 

reasonable mass. Any tangential motion would increase the 
inferred total LG mass as there would be a larger centrifugal 
force between the MW and M31 (which is not included in our 
analysis). This would tend to reduce model-predicted GRVs 
of LG galaxies, making it even more difficult to explain the 
observations. 

Including the LMC hardly affects our inference on Vc,q 
(Table 7). We used a prior on this parameter of 239zb5 
km/s (McMillan 2011). Our analysis slightly reduces its 
uncertainty. Based on the magnitude of this reduction, we 
conclude that our timing argument analysis independently 
constrains the speed of the LSR to within ~ 15 km/s. The 
combination of the prior with our work yields a best-htting 
LSR speed very close to that implied by the prior alone. 
This suggests that if we did not impose a prior constraint 
on Vc,Q, then we would hnd a most likely value for Vc,q that 
was within ~ x 0-5 ~ 5 km/s of 239 km/s. Thus, there 

is no tension between the LSR speed preferred by our anal¬ 
ysis and that preferred by McMillan (2011) on independent 
grounds, although our analysis is ~ 3 times less accurate in 
this regard. 

It is possible for massive satellites of M31 to cause a 
similar kinematic correction to its adopted HRV ^ . Our anal¬ 
ysis places a high statistical weight on M31 by requiring our 
models to match its HRV fairly well. However, our previous 
results were almost unaffected if we used the same value of 
O'extra for M31 as for other LG galaxies (the inferred value 
of oextra differed by ~ 2 km/s). This remains true if we 
include the kinematic effect of the LMC in our analysis. 
In other words, our results are not much different if we 
treat M31 in exactly the same way as other LG galaxies. 
Consequently, even a substantial alteration to the HRV of 
M31 should hardly affect our analysis as it involves 31 other 
galaxies. 

The case is different for the LMC because including 
a kinematic correction for it alters the predicted HRV of 
every galaxy in the LG rather than just the observed HRV 
of one galaxy. As the former does not much affect our overall 
conclusions, we suspect the same is true of the latter. 

4.5 Interactions With Massive Satellite Galaxies 

Models invoking gravitational slingshot encounters near the 
MW/M31 to hing out galaxies at high speed seem to share 
an important feature with the data: in Figure 8, most of the 
galaxies (20/32) are going outwards faster than predicted 


^ though see Peharrubia et al. (2016) for why such corrections 
are likely very small, even for the brightest M31 satellites 
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by the best-fitting model. This trend is perhaps clearer in 
Figure 9. 

We already included gravity from the MW and M31. 
One possibility not previously considered is that their satel¬ 
lites could have interacted with what are now non-satellite 
galaxies in the LG. For example, the Triangulum galaxy 
(M33) might be able to expand the region around M31 with 
a disturbed velocity field. This is possible via gravitational 
slingshot encounters with M33, using energy from its orbital 
motion around M31 to fling out material at high speed. 

Considering that M33 rotates at ~100 km/s (Corbelli 
2003), it can’t have affected the motion of a passing object by 
much more than this without merging with it. Thus, an im¬ 
portant issue with such a scenario is whether it can explain 
the fast outward motion of galaxies like DDO 125. Not only 
would it have to reach its present position several Mpc from 
M31, it would also have to possess sufficient kinetic energy 
to move at its present high velocity. Even if we neglect the 
retarding effect of gravity from the MW and M31, Hubble 
drag alone would mean that a peculiar velocity of 100 km/s 
today needed to have been 300 km/s at redshift 2, a plausible 
time for the interaction considering how far DDO 125 is from 
M31 (Figure 3). 

Moreover, one expects only a small fraction of the ma¬ 
terial in the LG to have interacted with Triangulum in the 
narrow range of impact parameters that lead to a large im¬ 
pulse but avoid a merger. Some of the material that was 
unaffected by M33 would no doubt have interacted with the 
LMC or with other massive satellites. Still, we find it hard 
to believe that such interactions would be as likely or as 
strong as required to fit the observations. Achieving both 
simultaneously does not seem feasible. 

Our models did not have particles starting too close to 
the MW or M31. We mapped the gravitational potential at 
t = ti (Equation 43) and assumed that all material below a 
certain level (i.e. with U < Uexc) had gone into one of these 
galaxies. Eor most parameters, this region did not split into 
separate regions around each galaxy but was a single region 
encompassing both. Test particles were not allowed to start 
within it. 

It is possible that pockets of material within this ‘ex¬ 
cluded region’ did not get accreted by the MW or M31. 
Starting closer to one of these galaxies, this material might 
be more likely to interact with one of their satellites. How¬ 
ever, it is unclear how such interactions could have been 
strong enough to explain the observations as the material 
would also have a deeper potential well to climb out of. 

4.6 Interactions with the MW and M31 

Other than the MW and M31, none of the objects in the 
LG seem heavy enough to impart a sufficiently large im¬ 
pulse on our target galaxies. However, our models already 
include gravitational slingshots caused by the MW & M31. 
Such encounters provide a way of extracting energy from the 
motion of these galaxies and putting it into the motion of a 
less massive third object. 

In principle, Andromeda can exert a large impulse on a 
passing object — perhaps up to twice Andromeda’s rotation 
speed. Therefore, it might be able to exert an impulse of as 
much as ~ 450 km/s on an object which approached closely 
enough yet avoided merging. Eor such an interaction to help 


explain the observations, the scattered object must have 
been fast-moving relative to M31 even when the two were 
far apart Otherwise, even fully reversing the small relative 
velocity ‘at infinity’ would only lead to a small impulse. 

In our simulations, the MW and M31 have never been 
moving very fast (Eigure 1). Their relative motion has usu¬ 
ally been slower than at present (~110 km/s, van der Marel 
et al. 2012a). It is difficult to achieve an impulse much ex¬ 
ceeding the motion of the massive body. Supposing the MW 
was moving at ~90 km/s, a small fraction of the material in 
the LG received an impulse of perhaps that much.^ 

The effect of Hubble drag then reduces the peculiar ve¬ 
locity gained in this way. So also does the gravity of M31 
(except for particles between the MW and M31, a region in 
which none of our target galaxies lie). Thus, the region in 
which the velocity field is disturbed by interactions with the 
MW/M31 only goes out to Mpc from the LG barycentre 
(Eigure 3). At higher LG masses, this region is somewhat 
larger: its linear size Mi^. It is difficult to see how this 
region can be made to encompass the whole LG. 

Prior to the start of our simulations, the MW—M31 
mutual gravity would not yet have had much time to retard 
their motion. Thus, we can assume that they were tracing 
the cosmic expansion, with mutual separation d{t) oc a{t). 
At these times, the Hubble parameter ^ ~ a “2 (Equation 
36) and so we expect the velocities of the galaxies to behave 
as a“ 2 . As a result, an interaction with a passing dwarf 
galaxy could lead to a maximum impulse on it that depends 
on the encounter time as ~ a (t) ~ 2 , This means that en¬ 
counters of the MW/M31 with LG dwarfs at very early times 
may have been very powerful. 

However, due to Hubble drag, the effect on the present 
peculiar velocity of the dwarf galaxy would be reduced by a 
factor of a at the time of the encounter.^ This means that 
very early encounters between the MW/M31 and LG dwarf 
galaxies should hardly affect our analysis, even if they were 
very strong. This is probably why our results changed very 
little when we altered the start time of our simulations to 
correspond to redshift 14 rather than 9 {cextra decreased by 
~ 1 km/s when using the earlier start time). 

It is unlikely that the MW and M31 existed at earlier 
times. This makes it difficult to argue that early encounters 
between LG dwarfs and the MW/M31 are responsible for 
the anomalously high HRVs of some distant LG galaxies. 

Therefore, one possible solution is to suppose that the 
MW and M31 were moving much faster than in our model 
a few Gyr after the Big Bang. We mentioned in Section 1 
that they might indeed have done so. In MOND, they would 
have undergone a close flyby ~9 Gyr ago (unlike the AGDM- 
based trajectories used in our models, Eigure 1). The relative 
speed at the time of closest approach would have been ~600 
km/s (Zhao et al. 2013). One could suppose that the MW 
was moving at 400 km/s and M31 at 200 km/s. Any passing 
dwarf galaxies would then have received a large impulse. 
Hubble drag would reduce peculiar velocities gained in this 
way by only a factor of ~3. Thus, the MW—M31 relative 

^ In a logarithmic potential, an extremely eccentric orbit has an 
angle of ~240° between apocentres, meaning that the deflection 
angle is ~60° rather than 180°. 

^ Figure 4 suggests that a factor of might be more accurate. 
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speed in this model seems about right to explain the motions 
of the LG galaxies with high AGRVs. 

Of course, one does not expect all of our target galaxies 
to have received an impulse quite this large. Only some of 
the material in the LG would have closely approached the 
MW or M31 at a time when their relative speed exceeded 
e.g. 400 km/s. This material might then get flung outwards 
and become an observed galaxy. It is also possible for the 
material to later merge with a galaxy that never strongly 
interacted with the MW or M31. Depending on the mass 
ratio of such a merger, it might leave the resulting object 
with a GRV only a little above that of the unperturbed 
galaxy before the merger. 

An interesting aspect of the observations which may 
point towards this scenario is apparent in Figure 10: the 
galaxies with the greatest excess radial velocity relative to 
AGDM all seem to be towards the edge of the LG. This may 
be because those objects which were flung out at higher 
speeds are now further away from the MW and M31. 

In this scenario, the high AGRV galaxies were all 
(roughly) at the same place at the same time: close to the 
LG barycentre when the MW—M31 flyby occurred. Thus, 
the magnitudes of the GRV discrepancies can be used to esti¬ 
mate the time of the flyby. The highest AGRV galaxy in our 
sample is DDO 99, with AGRV = 100 km/s and a distance 
of ~3 Mpc.^ Assuming objects at this distance would nearly 
follow a pure Hubble flow in AGDM (e.g. bottom panel of 
Figure 3), its radial velocity should he H^d 200 km/s with 
respect to the LG bary centre. Neglecting projection effects,^ 
the actual radial velocity is ~50% larger. This corresponds 
to an elapsed time since the flyby of ~ | the age of the 
Universe, i.e. ~9 Gyr ago. 

Interestingly, this is also when there appears to have 
been a sudden perturbation to the disc of the MW which cre¬ 
ated its thick disc (Quillen V Garnett 2001). There is some 
circumstantial evidence that this perturbation was tidal in 
nature rather than a process internal to the MW (Banik 
2014). Although a minor merger is an obvious possibility, 
it would leave accreted stars. As the accreted galaxy must 
have been reasonably massive compared to the MW to create 
its thick disc, these accreted stars should have characteristic 
properties. Recent attempts to find such stars have not found 
any (Ruchti et al. 2015). This might be an indication that 
the thick disc was created by a close flyby of another massive 
galaxy ~9 Gyr ago rather than a minor merger at that time. 

In this respect, it is interesting to estimate when a 
MW—M31 interaction might have occurred if MOND were 
the correct description of nature. Applying this theory, it 
can be shown that the time from apocentre to pericentre is 
given by (equations 15 and 29—30 of Zhao et al. 2010) 


At 



dapo 

^OMGMa^ 




1 - 1/2 


dx 


where 


( 66 ) 


Note that the relevant mass M is the total baryonic 
mass of the LG. The non-linear nature of the theory reduces 
the force between two particles with the same total mass if it 
is distributed more equally. Even assuming (conservatively) 
equal MW and M31 masses (leading to the factor of 0.61) 
and using a very low estimate for M of 10^^M©, we would 
get At ~ 7.2 Gyr.^ The MW and M31 are slightly past their 
apocentre now, but it is still clear that they must have had 
a past close flyby in this theory. 

Equation 66 neglects several complications which arise 
in MOND. Most important is the external field effect 
(Bekenstein V Milgrom 1984; Milgrom 1986). This arises 
because MOND is an acceleration-dependent theory. As a 
result, a constant external gravitational field acting upon a 
system weakens the self-gravity of objects within the system. 
This effect is approximately taken into account in the work 
of Zhao et al. (2013). 

We are currently undertaking more accurate timing ar¬ 
gument calculations in the context of MOND. Preliminary 
results indicate that the LG mass implied by the MOND tim¬ 
ing argument is consistent with baryons only, if one requires 
a past close approach between the galaxies. This flyby needs 
to have been ~8—9 Gyr ago. It is not feasible to construct 
trajectories of the MW and M31 in MOND which avoid such 
a close approach or have ^2 such events. 

In MOND, the longer range nature of gravity means 
that we need a more careful treatment of objects outside the 
LG (Table 1). It is possible that some of the anomalously 
high outwards velocities found in this work are due to LG 
galaxies — especially those close to perturbers — falling in 
towards them. Indeed, including tides from Gentaurus A im¬ 
proved the fit to observations somewhat, even in the context 
of AGDM (Eigure 14). This effect might be further enhanced 
in modified gravity scenarios. 

Interestingly, the discrepancy does seem to be higher 
for target galaxies closer to M81 or to IG 342 (Eigure 15). 
We argued that the apparent correlation could not be due 
to gravity from these perturbers as this would imply that 
they affected velocities of target galaxies more than is rea¬ 
sonable. However, this argument likely breaks down under 
a different law of gravity. In this case, tides might well be 
more significant than we assumed. 


4.7 Comparison with Peharrubia et al. (2014) 

Our results are broadly similar to the recent study con¬ 
ducted by Peharrubia et al. (2014), hereafter P14. Those 
authors also favour a low and a similar LG mass. We 
found a slightly higher value of a extra (45lg km/s instead 
of 35^4 km/s), though the estimates are marginally in agree¬ 
ment. 

However, our investigation greatly strengthens the con¬ 
clusions reached by P14. We used an axisymmetric model for 
the Local Group rather than a spherically symmetric one. 
As our target galaxies are often not much further away than 


^ As this is the highest AGRV out of 32 galaxies, the true value 
is likely smaller than the 110 km/s obtained using nominal values. 
^ on the sky, DDO 99 and M31 are almost at right angles (Figure 
3), so our conclusions should not be much affected by uncertainty 
in the motion of the MW due to uncertainty in 


^ We used cosmological parameters as in Table 2 and an apoc¬ 
entre distance dapo = 1 Mpc, slightly less than in ACDM due to 
the stronger gravity. For o-q, we used 1.2 x 10“^° m/s^ (McGaugh 
2011 ). 
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the MW—M31 separation, a spherically symmetric gravi¬ 
tational field may be a poor approximation. For example, 
gravitational slingshot encounters with the MW/M31 rely 
on a time-dependent non-spherical potential. Without mod¬ 
elling either of these effects, it would be difficult to draw 
reliable conclusions about whether these close encounters 
might have left an imprint on the present motions of target 
galaxies. 

Even if one could be sure that such encounters were not 
important, a two-centred potential has other subtle conse¬ 
quences. A trajectory which initially went orthogonal to the 
MW—M31 line from the point halfway between them; would 
curve towards the heavier galaxy (almost certainly M31). 
This would tend to increase the HRV of the target galaxy 
(e.g. bottom panel of Figure 5). Curvature of test particle 
trajectories seems to be important even at quite large dis¬ 
tances from the MW and M31 (top panel of Figure 3). The 
process is more significant if the MW and M31 masses are 
very unequal, which definitely seems to be the case (Figure 

7 ). 

By using the same list of target galaxies as P14, we avoid 
targets too close to any of the major perturbers relevant to 
our analysis (Table 1). However, tides from these objects 
must affect our results at some level. We directly include 
the most massive perturber (Centaurus A), exploiting its 
location almost along the MW—M31 line (Section 2.2.2). 
We think this greatly improves our model. We conduct a 
more thorough investigation of tides raised by other objects 
(Section 4.3) and consider the effect of the Large Magellanic 
Cloud in some detail (Section 4.4). 

Our initial conditions are handled differently to P14. We 
use cosmological initial conditions (Equation 37) because of 
observations indicating very low peculiar velocities at early 
times (Planck Collaboration 2015). P14 used a procedure 
involving non-cosmological initial conditions which does not 
seem entirely physical (see their section 3). 

We added an extra term to our equation of motion 
(Equation 40) to account for cosmological expansion. The 
idea is to recover r (x a (t) in the absence of inhomogeneities. 
A similar approach was used by P14. However, they did 
not include the deceleration to the expansion rate caused 
by matter, leaving only the acceleration caused by dark 
energy. This implicitly assumes that the rest of the LG is 
empty. As we did not make this assumption, we suspect that 
the predicted HRVs in our investigation are lower. With a 
present dark energy fraction of ~ 0.7, we expect a difference 
of ~ (l — a/0?7) HqcI for a target at distance d. Assuming 
a typical distance of 2 Mpc, this suggests a difference of 
~20 km/s. Because observed HRVs tend to exceed predicted 
ones, it is unsurprising that our estimate of a extra is higher. 

In Section 4.1, we discuss how much mass might actually 
be in the RLG. Here, we show that, although it is possible 
to have an empty RLG, this is an extreme case. It is also 
possible for it to have even more mass than we assumed. 
Our assumption is an intermediate case. Nonetheless, we 
performed calculations assuming an empty RLG and showed 
that this reduces a extra by ~7 km/s. 

We argued that M31 should be treated specially in that 
one should use a lower value of a extra for it than for other LG 
galaxies. This forces our models to match the HRV of M31 
very well, restricting which models are viable and thus rais¬ 
ing (Textra ‘ The cffcct is Substantial for our analysis without 


Gen A: the most likely value of a extra is raised from 46 to 54 
km/s. However, in our more realistic analysis including Gen 
A, a extra is Only affected by ~2 km/s. Thus, although we 
recommend treating M31 specially due to its much higher 
mass than LG dwarf galaxies, our overall conclusions are 
little altered if one does not do so. The parameter most 
affected seems to be Mi, which is lower if M31 is not treated 
specially. This is also apparent in Figure 13 of P14, especially 
if one imposes an independent constraint on the LSR speed 
(e.g. McMillan 2011). 


5 CONCLUSIONS 


We performed a careful dynamical analysis of the Local 
Group (LG) to try and explain the observed positions and 
velocities of galaxies within it. The LG was treated as a 
collection of test particles and two massive ones — the Milky 
Way (MW) and Andromeda (M31) — which we put on a 
radial orbit. We added a third massive particle to repre¬ 
sent Gentaurus A, which is very close to the MW—M31 line 
(Table 1). All particles were started moving outwards from 
the centre of mass of the LG with speeds proportional to 
distance from there. Thus, they all started on a pure Hubble 
flow with no peculiar velocity (Equation 37). 

A wide range of possible masses for the MW and M31 
was investigated using a grid method (Table 2). Each time, 
we got the final MW—M31 and MW—Gen A distances to 
match their observed values. We also got a test particle 
trajectory to end at the same location as each observed LG 
galaxy. This gave a model-predicted velocity, whose line-of- 
sight component (the HRV) was compared with observa¬ 
tions. 

The best-fitting total LG mass is 4.33lo'32 x lO^^M©, 
with 0.14 ± 0.07 of this accounted for by the MW and the 
rest by M31. There is almost no tension between the Lo¬ 
cal Standard of Rest rotation speed estimated by McMillan 
(2011) and our analysis (Table 7). 

However, even in the best-fitting model, there was a 
poor match between observed and model-predicted HRVs. 
Thus, we tried to quantify the extra astrophysical noise 
O'extra that the observations imply. To do this, we added it 
in quadrature with the other known sources of error, which 
are all observational. Oextra can be constrained using 


P{Mode\\HRVohs) oc -e 

<7 


r HRV^u,-HRV^, 


(67) 


Raising oextra — and thus a — will eventually cause the 
probability to decrease^. In this way, we found that oextra — 
45 . 115'7 km/s (Figure 7). This rather high value is partly 
due to some galaxies receding from the LG even faster than 
a pure Hubble flow (Figure 11). This is despite the attractive 
gravity of the MW and M31. 

We expect interactions between LG dwarf galaxies to 
have contributed somewhat to oextra, but only at the ~15 
km/s level. This is because they have fairly low rotation 
velocities/internal velocity dispersions (Kirby et al. 2014), 
limiting how much they can influence each other gravita¬ 
tionally. If we consider just those galaxies with HRVs below 
the predictions of the best-fitting model (blue area in Figure 


1 when cr > \HRVobs - HRVnoded 
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9), then we see that they can indeed be described quite well 
by a Gaussian of this width. 

One might expect the same to be true for LG galaxies 
with HRVs that exceed the model predictions. However, this 
is not true (red area in Figure 9). Thus, there is a systematic 
trend for radial velocities to be higher than we expect. 

We considered tides from objects outside the LG at 
some length (Section 4.3). The most relevant perturbers are 
given in Table 1 . A correlation is apparent whereby the most 
problematic galaxies are closest to these perturbers (Figure 
15). Thus, we constructed a simplified model to estimate 
how much they might have affected the GRVs of the most 
discrepant galaxies. Due to a combination of projection ef¬ 
fects and large distances from the perturbers (> 2 Mpc), 
we consider it unlikely that tides from IG 342 and M81 can 
reconcile our model with observations (Table 6). In fact, 
they would likely exacerbate the tension in several cases. 
This also seems to be true of tides raised by the Great 
Attractor (Figure 16): including these raises a extra by ~5 
km/s (Section 4.3.3). 

A local under-density may help to explain the observa¬ 
tions. We consider this possibility in Section 4.1 and show 
that this can reduce a extra by at most ~7 km/s. Increas¬ 
ing the Hubble constant Hq has a similar effect, which we 
consider in Section 4.2. This reduces a extra by km/s 
(Figure 14). Although is known quite accurately (Planck 
Gollaboration 2015), a higher value may be appropriate for 
the LG if there is a local under-density. 

Satellite galaxies of the MW can affect our analysis 
kinematically. In this regard, we considered the Large Mag¬ 
ellanic Gloud (LMG) in some detail (Section 4.5). Including 
the LMG can reduce a extra by at most 5 km/s (Figure 17). 
Thus, incorporating it into our models does not reconcile 
them with observations. 

Interactions of LG dwarf galaxies with massive 
MW/M31 satellites (e.g. M33) would probably be too weak 
to explain distant non-satellite galaxies moving outwards 
even faster than a pure Hubble flow. As for encounters with 
the MW and M31 themselves, this would naturally explain 
why LG galaxies tend to move outwards much faster than 
expected. However, this process is already included in our 
simulations — it seems to be too weak. 

It might have been more powerful in reality if the MW 
and M31 had been moving much faster than in our mod¬ 
els (Section 4.6). Although this would be very unusual in 
the context of AGDM, such fast motions arise naturally in 
certain modified gravity theories (Llinares et al. 2009). For 
example. Modified Newtonian Dynamics (MOND, Milgrom 
1983) leads to a past close encounter between the MW and 
M31 and a maximum relative speed of ~600 km/s (Zhao 
et al. 2013). 

One aspect of our results which may point towards this 
scenario lies in the distances to the objects with the highest 
radial velocities relative to the best-fitting AGDM model. 
These distances all seem to be close to the upper limit of the 
distance range probed by our sample (Figure 10). This may 
be because those objects which received the strongest gravi¬ 
tational boost from a close encounter with the MW/M31 are 
currently furthest away from the location of the encounter. 
Gonsidering the distances and velocities of such galaxies sug¬ 
gests that the encounter was ~9 Gyr ago. This is around the 
same age as the thick disc of the MW (Quillen & Garnett 


2001). It is also roughly when a MOND timing argument 
calculation suggests that a MW—M31 encounter took place 
(Zhao et al. 2013). 

We think it likely that the observed positions and veloc¬ 
ities of LG galaxies would be easier to explain if there was a 
past close MW—M31 flyby. This is possible only within the 
framework of a modified gravity theory. More work will be 
required to test such a scenario. 

The authors wish to thank Jorge Penarubbia for helpful 
discussions and the referee for their patience while we im¬ 
plemented their suggestions, which greatly strengthened this 
contribution. IB is supported by a Science and Technology 
Facilities Gouncil studentship. The algorithms were set up 
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